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ABSTRACT 


This  thesis  presents  a  model  which  simulates  the  scattering  from  a  fluid 
loaded  I-beam  and  the  resultant  behavior  due  to  fluid-structure  interaction. 
Chapter  I  gives  an  overview  of  the  problem  and  describes  the  characteristics 
of  the  solid  and  fluid,  the  aspects  of  periodicity,  boundary  conditions  and  the 
coupling  of  the  two  media. 

The  governing  equations  of  motion  are  scaled  in  Chapter  II.  In  Chapter 
III,  the  finite-difference  formulae  fcr  these  equations  are  derived,  as  is  the 
non-local  radiation  boundary  condition.  Difference  formulas  for  typical 
boundary  points  of  the  solid  and  corner  nodes  are  also  derived.  All  finite 
difference  formulae  used  are  presented  in  Appendix  C.  Chapter  IV  contains 
numerical  results.  Conclusions  are  drawn  and  areas  of  the  problem  that 
would  require  further  study  are  in  Chapter  V. 
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I.  PROBLEM  OVERVIEW  AND  DOMAIN  DESCRIPTION 


The  problem  we  consider  is  where  an  acoustic  pressure  wave  emitted  by 
an  active  sonar  impinges  on  a  submarine.  In  our  case  this  is  a  double-hulled 
submarine.  Active  sonar  works  on  the  principle  of  detecting  the  reflection  off 
a  solid  object  of  an  acoustic  pressure  wave  emitted  by  a  source,  by  which  one 
can  calculate  distance  and  direction  to  that  object  as  in  Figure  1,  but  instead  of 
concentrating  on  the  reflected  wave  we  look  at  the  interaction  between  the 
structure  and  the  incoming  wave.  This  generates  scattered  pressure  waves. 
The  scattered  pressure  waves  include  waves  which  decay  as  they  travel 
through  the  fluid  (evanescent  modes),  and  waves  which  do  not  (propagating 
modes).  Since  propagating  modes  do  not  decay  they  can  be  detected.  The 
main  thrust  of  this  thesis  is  to  determine  the  characteristics  of  the  propagating 
modes,  such  as  amplitude  and  energy.  The  optimum  situation  wouid  be  to 
perturb  the  double-hulled  structure  at  a  resonance  frequency  which  would 
increase  the  amplitude  of  the  propagating  modes  making  them  easier  to 
detect.  We  investigate  the  steady  state  behavior  of  the  propagating  modes  and 
the  resultant  shear  strain  field  in  the  I-beam.  At  steady  state  the 
characteristics  of  the  propagating  modes  stabilize  which  might  be  used  as  an 
acoustic  signature  of  the  structure.  In  addition  when  the  shear  strain  field 
reaches  steady  state  its  value  throughout  the  solid  may  be  used  to  isolate  areas 
of  the  I-beam  where  large  stresses  occur.  This  information  could  be  useful  in 
the  design  of  double-hulled  vessels. 
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To  reduce  the  model  to  one  where  we  can  apply  numerical  techniques  to 
simulate  its  behavior  we  must  make  simplifying  assumptions  which  are 
discussed  below. 


A.  ASSUMPTIONS 

•  The  source  of  the  pressure  waves  is  far  enough  away  so  that  the 
impinging  wave  can  be  approximated  by  a  normally  incident  plane 
wave. 

•  There  are  no  other  sources  present  such  as  might  occur  with  bottom 
bounce,  surface  reflection  and  the  like. 

•  The  area  of  interest-— AOI  as  shown  in  Figure  2  where  the  wave 
impinges  is  where  the  inner  and  outer  hulls  are  joined  or  connected  by 
a  supporting  spar  forming  an  I-beam  shaped  domain. 

•  The  I  beam  shaped  domain  is  a  uniformly  continuous  linear  isotropic 
elastic  medium  with  no  cracks,  welds  or  other  deformities. 

•  The  dimensions  of  our  AOI  are  such  that  any  curvature  of  the  surfaces 
can  be  ignored. 

•  The  center  spar  or  support  beam  occurs  at  regular  intervals  through  the 
structure  as  shown  in  Figure  3  allowing  us  to  truncate  the  domain  to 
the  left  and  right  of  center  spar  using  periodic  boundary  conditions. 

•  The  incident  wave  does  not  displace  the  submarine. 

•  The  cavities  A  and  B  in  Figure  2  enclosed  by  the  inner  and  outer  hulls, 
and  the  center  spar  is  void  and  contain  no  sources. 

•  We  are  only  interested  in  displacement  in  the  x  and  y  directions  as 
given  in  Figure  2  which  reduces  our  problem  to  two  dimensions. 

•  The  fluid  is  seawater,  the  solid  is  steel. 

Our  model  is  now  reduced  to  one  where  we  have  two  coupled  media — fluid 
and  solid,  and  the  characteristics  of  each  will  now  be  discussed  in  turn. 
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Figure  2.  Area  of  Investigation 


A.  THE  SOLID 


A  cross-section  of  the  double  hull  of  length  8a  is  shown  in  Figure  3.  Note 
a  is  a  scaling  constant  for  distance.  The  domain  is  essentially  infinite  in 
extent  in  the  x  (length)  and  z  (depth)  directions.  Since  it  is  of  uniform  shape 
in  the  z  direction,  we  can  neglect  this  coordinate  and  reduce  our  problem  to 
two  dimensions.  Our  area  of  investigation  is  outlined  in  Figure  3,  and  we  are 
now  faced  with  the  problem  of  providing  boundary  conditions  in  the  x 
domain. 


Figure  3.  Cross-section  of  Submarine  Hull 


We  have  a  normally  incident  plane  wave  impinging  over  one  boundary 
of  the  domain.  Thus  the  pressure  is  the  same  at  all  points  (there  is  no  phase 
shift)  along  the  fluid /structure  boundary.  Due  to  the  periodic  nature  of  the 


structure  we  expect  the  solid  to  behave  the  same  over  given  intervals  of  2a, 
that  is 

f(x,y,t)  =  f(x  +  2aN,y,t)  (1) 

where  N  is  an  integer  and  /  a  function  describing  the  state  of  the  solid  such  as 
displacement,  velocity  and  so  on.  Using  this  periodicity  we  can  truncate  our 
domain  at  a  and  -a  and  use  the  following  periodic  boundary  conditions 

f{-a,y,t)  =  f(aiy,t)  (2) 


and 


/*(-a,t/,f)  =  /*(a,y,f) 


(3) 


where  k  is  a  positive  integer  and  can  signify  the  derivative  with  respect  to  x,  y 
or  i,  depending  on  the  function  /. 

The  solid  has  now  been  reduced  to  a  two-dimensional  linear  isotropic 
elastic  medium  whose  governing  equation  of  motion  for  points  interior  to 
the  domain  is  given  by  the  plane  strain  elastic  wave  equation  which  is 


f  d2u  B2u) 


a**V 


+  (A+/i) 


!i2u  3V 


Bx 


2  +  dxdy 


B2u 

~Psdt2 


(4) 


f  d2v 

dx2  By2 


^  B?'v  t  d2u 
Bx 2  BxBy 


+  (A  +/i) 


B2v 

~PsW 


(5) 


u  and  v  are  displacements  in  the  lateral  (x)  and  transverse  (y)  direction,  p  and 
X  are  Lam6  constants,  and  ps  is  the  density  of  the  solid.  The  I-beam  is 
deformed  due  to  the  imposed  fluid  pressure.  Balancing  normal  forces  at  the 
fluid/ solid  interface  we  obtain  the  boundary  condition 
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_  _  ildu  ,  dv),  rj'dv  _  total 


(6) 


where  ptotal  is  the  total  pressure,  and  Xyy  the  normal  stress  component.  All 
other  boundary  conditions  are  either  traction  free  or  periodic  and  are 
summarized  below. 


Ux  =  01 


T^=°, 


Xyy  -  01 

T^  =  0 


[  on  surface  DH  and  El  in  Figure  4, 


(7) 


on  CD,  EF,  GH,  IJ,  and  KL  in  Figure  4, 


(8) 


r  =  -ntotal 

Tyy  P 


A  D  iw  /vt  i  ««r\  A 

f  vu  nL/  iii  i  -x/  uni*. 


fc\\ 

\y) 


periodic  ] 
conditions 


elsewhere. 


(10) 


%xx,  tyy  and  %y  are  the  normal  and  shear  components  of  stress  and  are 
represented  by 


‘yy 


T XX  -  ^ 


f  9u  91’^ 

'dv' 

1&+¥J+^ 

9u  dv  |  _  f9u\ 

*+*  J+n&.l'“d 


9u  9^ 

,9y H  9* 


(11) 


(12.) 


(13) 


7 


Figure  4  is  the  unsealed  domain,  in  which  s,  and  k  are  constants  used  to  vary 
the  size  of  the  cavities  A  and  B  of  Figure  2.  Note:  U  <  s  <  1,  0  <k<  1. 


Figure  4.  Unsealed  Domain 


B.  THE  FLUID 

We  are  interested  in  the  scattered  pressure  waves  generated  by 
perturbations  at  the  fluid/solid  interface.  The  partial  differential  equation 
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governing  the  propagation  of  these  waves  in  the  fluid  is  the  linear  two- 
dimensional  wave  equation  given  by 


1  d2p  _  d^p  d2p 
’cjHt2~l)x2+  dy2 


Q4) 


where  Cf  is  the  speed  of  sound  in  the  fluid. 

These  perturbations  arise  when  a  normal  plane  wave  of  magnitude  P  and 
frequency  (0  impinges  on  the  surface  of  the  solid,  causing  it  to  deform.  Since 
the  pressure  wave  is  of  normal  incidence  (there  is  no  phase  shift  at  points 
along  the  surface  of  the  solid)  and  the  solid  is  a  periodic  structure  we  can 
expect  the  behavior  of  tne  scattered  pressure  waves  to  be  the  same  in  given 
intervals  of  2a,  thus  we  can  use  periodic  boundary  conditions  in  the  same 
manner  as  was  used  for  the  solid. 

Ideally,  if  the  solid  were  acoustically  hard  our  total  pressure  would  only 
have  two  components,  incident  and  reflected.  In  reality  the  solid  is  perturbed 
by  the  incident  wave,  generating  scattered  pressure  waves  which  propagate 
out  into  the  fluid  domain.  The  total  pressure  in  the  fluid  can  be  represented 

ty 

ptotal  =  ^incident  +  ^reflected  +  ^scattered  (14) 

where  ^incident  =  ^reflected  =  e'kP'l[  and  ^scattered  js  to  be 

determined.  To  avoid  cavitation  at  the  fluid-solid  boundary  we  employ  the 
inviscid  form  of  the  Navier-Stokes  equation  which  we  refer  to  as  the 
compatibility  condition  and  is  given  by 
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(15) 


Vj  _ 

dy  !y=o  ^  Bt 2 


where  pf  is  the  density  of  the  fluid.  Note  that  only  the  scattered  term  of  the 
pressure  is  included  in  this  equation  since 


Bp1 


dy 


=  -ikfe 


-it 


R 


and 


y=G 


dy 


y~o 


-  ikfe 


(16) 


thus 


=0. 

,  dy  By  }y= o 


(17) 


The  scattered  pressure  waves  generated  at  the  fluid/solid  interface  are 
composed  of  propagating  (non-decaying)  and  evanescent  (decaying)  modes 
which  must  be  allowed  to  propagate  off  to  infinity  in  the  positive  y  direction. 
To  do  this  we  must  employ  a  non-local  radiation  boundary  condition  whose 
implementation  will  be  discussed  in  greater  detail  in  the  discretization 
section  of  this  paper. 
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H.  NON-DIMENSIONALIZATION  OF  THE  GOVERNING  EQUATIONS 

OF  MOTION  AND  BOUNDARY  CONDITIONS 


The  problem  addressed  models  the  propagation  of  a  scattered  pressure 
wave  in  two  dimensions.  This  is  described  by  the  two-dimensional  scalar 
wave  equation 


(18) 


where  p  is  the  pressure  and  Cf  the  speed  of  sound  in  the  fluid.  To  facilitate 
implementation  and  to  free  ourselves  from  the  requirement  of  using  a  given 
system  of  measurement  such  as  metric  or  imperial  we  scale  or  non- 
dimensionalize  Equation  18  as  follows.  Let 


1  d2v  __  9 2p  d2p 
cj  dt2  dx2  dv2 


(19) 


represent  the  unsealed  wave  equation.  We  now  use  the  following 
relationships: 

t  =  cot;  x  =  ~;  y--;  p  =  ~.  (20) 

a  a  P 

Here  (o  is  the  scaling  constant  for  time  and  the  frequency  of  the  incident  plane 
wave,  a  is  the  half  length  of  the  I-beam  and  the  scaling  constant  for  distance 
and  P  is  the  scaling  constant  for  pressure.  Note  that  when  taking  derivatives 
with  respect  to  x  we  get 
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and 


(21) 


A  =  IA. 

dx  a  dx 


dx2 


1  dl 
a2  dx2 


This  also  holds  true  for  derivatives  with  respect  to  y,  and  a  similar  relation 
results  when  taking  derivatives  with  respect  to  t.  With  the  above 
relationships,  Equation  19  is  now  written  as 


afP  d2p  _  P  <rp  P  d  2p 
1f'dtT~'aId^  +  ~a2dy2' 


(22) 


Cancelling  common  factors  and  multiplying  Equation  22  by  a 2  reduces  it  to 


ora2  d2p  _  d2v  crp 


c2  dt2 


*x2  +  By2  ‘ 


(23) 


Defining 


,  l oa 

kf  =  T,' 


Equation  23  is  now  written  as 


J a*Lj2JL+iiL. 

*  dt2  dx2  dy2 


(24) 


(25) 


The  elastic  wave  equation,  which  is  a  vector  wave  equation  is  given  by 

d2u 


0 d2u  d2u^ 

vi?  +  ^2' 


+  (A  +/i) 


1  d2U  [  d2V  ^ 

dx2  dxdy 


=  Ps 


at" 


(26) 


fd2  v  d2v) 

- 1 - 

a^2  ^2 


+(*+/*) 


(  d2v 


d2u  ^ 


(ay2  a*ay 


d2V 

dp 


~Ps\  2  ■ 


(27) 


The  following  relationships  allow  us  to  write  Equations  26  and  27  in  a  more 
convenient  form. 


—  =  Cj  and 
Ps 


X  +  2fi 
Ps 


(28) 


The  constants  cl  and  cj  are  the  lateral  and  transverse  velocities  of  the  solid. 
The  unsealed  forms  of  Equations  26  and  27  are  now  written  as 


CT 


^ d2U  d2U ^ 
ydx2  +  df 


( d2u  a2-  ^ 


'V 


dxdy 


d2u 

W 


(29) 


0 

Cf 


^  d2v  d2o^ 
+ 


-2 


dx 


dy2 


(el-4) 


(  d2v 


■  4*  — 


d2u  ^ 


dy 2  dxdy 


d2v 

W' 


(30) 


Wc  use  the  same  scaling  relationships  as  before  for  x,  y  and  t  (Equation  20)  as 
well  as 


u 


and 


(31) 


which  are  the  scaling  relationships  for  the  displacement,  and  rewrite 
Equations  29  and  30  as 


2 ( D  d2u  D  d2u)  12  2\(  D  d2u  D  d2V  „  2^u 


4 


(  r-»  ^2, 


D  d^v  D  /  2  2/  D  d2v  D  dlu 

^a2  dx2  a2  dy2  J  3  \a2  dy2  a2  dxdy 


2.,  ^ 


n  2^ 
=  Deo  — =-• 


dP 


(33) 


Defining 


,  coa  ,  ,  cm 

kL=—  and  kT=~r 


(34) 
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and  dividing  Equations  32  and  33  through  by  o2  and  cancelling  common 
factors  they  reduce  to 


1  f  d2u  d2u 


i.2  3^2^" 


.2  +  l2 


kf{dx2  dy1 )  {kf  kf  j(dx2  3x3y J  dt 


1  1  Y  d^u  d2v  )  d2u 


1  f  d2ZJ  32d1  f_l _ 1  Y<A>  32u  )  rj2V 

if  la? + V  J + liF  ’  if  w +  r  s'7' 


Collecting  like  terms  gives  us  the  final  form 


1  d2u  1  d2u  1  _  1  Y  d2v  _  d2u 
kf  dx2  kj  dy2  [kf  kfj^xdy J  dt2 

i  d2v  i  d2v  (  \  iY32«1  az» 

kf  dx2  kfdy2  *r,  dxdy ^  dt2 


The  surfaces  in  contact  with  free  space  are  stress  free.  (The  fluid /solid 
boundary  is  dealt  with  separately).  Therefore  the  stresses  xxx  and  xXy  on 
surfaces  El  and  DH  and  %  and  xXy  on  surfaces  CD,  EF,  GH,  IJ  and  KL  of 
Figure  5  are  zero. 

These  components  of  the  stress  tensor  can  be  written 

'du  dv}  n  .OQX 

+  =0  (39) 

^  I9y  dxl 


,  f  du  dv  |  _  du 


Jar  an  ,  as  „ 

^-WayKaT0 
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where  x,  y~,  u  and  v  are  the  unsealed  components  of  distance  and 
displacement  and  fi  and  A  are  the  Lain£  constants.  We  will  now  scale 
Equations  39  through  41  in  turn. 

For  Equation  39  using  the  same  scaling  relationships  as  before  (see 
Equations  20  and  31)  gives 


15 


(Ddu  Ddv)  A 

P  ~ T-  +  —  3“  =  °- 
^  a  dy  a  dxj 


Cancelling  common  factors  reduces  it  to 

du  dv 

%  =  ~dj  +  dx  =  °'  (43) 

For  Equation  40  we  first  divide  through  by  ps  the  density  of  the  solid,  and 
using  the  same  scaling  relationships  we  get 

X  (Ddu  Ddv)  2fi  D  du  _ 
ps{a  ox  a  dy)  ps  a  dx 


We  know  that 


Similarly  it  can  be  shown  that 


p! Ps  —  cj- 


X  2  2 

=  cL-  2cj. 


Using  Equations  45  and  46,  substituting  into  Equation  44,  and  cancelling 
common  factors  gives 


£i._2Y|l+|£)+2^  =  o. 


cj  hdx  ty)  dx 


Using  the  previous  definitions  ot  ki  and  kj  it  can  be  shown  that 
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(49) 


£l  =  *r 

CT  kl 


and  when  used  in  Equation  48  it  reduces  to  our  final  form 


kl  du 


du 


Following  the  same  procedure  for  Equation  41  its  final  form  is 

'4-2l|u4  |^=o. 

Ui  )ix  kl  % 


(50) 


(51) 


At  the  interface  the  fluid  cannot  exert  a  force  tangential  to  the  boundary, 
hence  the  shear  component  of  stress  is  zero  there,  and  is  given  by 


'•xy 


f  "\T7  \ 

■  uu  uu  * 


=  0. 


(52) 


For  the  normal  component  of  stress  at  the  interface  we  refer  back  to  Equation 
6,  Section  B  and  Equation  15  in  Section  C  of  Chapter  I  to  see  that  the  normal 
component  of  stress  is  given  by 

.(du  dv')  .  du  t  i  r?  c \ 

’"■xl»  +  asJ+:v*‘"vr+r+'r>  (53) 


where  the  superscripts  I,  R  and  S  signify  incident,  reflected  and  scattered 
pressures.  The  incident  and  reflected  pressures  are  given  by 

p‘  -'■*!*■*  and  pR=a^-" 

Txy  is  scaled  the  same  as  before.  To  scale  Tyy  properly  we  will  need  the 
compatibility  condition,  whose  unsealed  form  is 
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(54) 


.§e! 

3y 


y  =  0 


d2v 


We  refer  back  to  Equation  17,  Chapter  I,  Section  C  to  see  that  only  the 
scattered  pressure  need  be  considered  here. 

Scaling  Equation  54  we  get 


Pdps 


a  3 y 


y  =  0 


=  pf(o 


D 


82n 

0f2' 


(55) 


(56) 


tquation  06  gives  us  a  convenient  choice  for  the  scaling  constant  for  pressure 
of 


P  =  oraDpf 


(57) 


and  when  substituted  back  into  Equation  56  cancels  as  a  common  factor  to 
give 


dp*  _d2v 
3y  y  =  0  dt2' 


(58) 


Dividing  Equation  53  by  ps  and  scaling  gives 


AD 

Psa 


rdu  ,  dv)  .  2pD  dv  _  -P/„/  ,  „r  ,  „s\ 

+  *  +P  +PL=i 


ps«  dy  ps 


y=0 


(59) 


Substituting  the  value  of  P  from  Equation  57  we  obtain 
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(60) 


A  D  du  Bv 
psa{dx  By 


2/iD  3t> 
p$a  By 


-5«!44p! 

Ps 


or 


f  3»  +  3:;^ 

+  2^  = 

-ekjlp1  +pR  +ps) 

U  J 

\Bx  By j 

'  v  r  r  / 

y  =  0 


where  e  =  p//ps. 

Evaluating  the  incident  and  reflected  pressures  at  y 
pi  =pR  =  e-P,  and  when  substituted  into  Equation  61  gives 


'*£_  2 


( Bu  Bv ^ 


+  2— =  -2 ek£e~u  -  ek$ps . 

<>y 


(61) 


0  we  get 


(62) 
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III.  DERIVATION  OF  THE  FINITE  DIFFERENCE  FORMULAE  FOR  THE 

GOVERNING  EQUATIONS 


Throughout  the  derivation  of  the  finite  difference  approximations  we 
identify  gridpoints  of  the  scaled  domain  with  the  subscripts  i  and  j  and 
superscript  n.  Variation  in  the  x  direction  is  denoted  by  i,  1  <  i  <  N,  where  N 
is  the  total  number  of  subdivisions  (since  we  are  using  a  square  grid  the 
number  of  subdivisions  in  the  y  direction  is  the  same),  variation  in  y  with  j,  1 
<  j  <  N,  and  time  with  n,  0  <  n  <  <».  Lower  case  indices  denote  varying 
quantities,  while  upper  case  letters  denote  fixed  quantities.  would  be  the 
value  of  the  scattered  pressure  for  all  values  of  i  at  the  grid  level  y  =  KAy,  thus 
is  a  vector,  while  is  an  element. 

As  '  discussed  above  we  use  the  letter  i  to  denote  variation  in  the  x 
directs  .  his  is  a  common  convention  and  we  do  not  wish  to  deviate  from 
it.  To  avoid  confusion  with  the  complex  quantity,  ^-1,  also  commonly 
denoted  by  i,  we  state  the  following  rule,  that  whenever  the  letter  i  appears  as 
a  superscript  it  denotes  the  complex  quantity  y[~l,  and  when  appearing  as  a 
subscript  it  is  an  index  denoting  variation  in  the  x  direction. 

A.  FINITE  DIFFERENCE  APPROXIMATIONS  FOR  THE  EQUATIONS 

GOVERNING  THE  BEHAVIOR  OF  THE  FLUID 

The  scaled  domain  of  consideration  as  shown  in  figure  6  has  periodic 
boundary  conditions  applied  at  x  =  1  and  -1  and  a  radiation  boundary 
condition  for  the  propagating  modes  at  y  =  2.  To  model  the  two-dimensional 
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Periodic 

Boundary 

Conditions 


P(W) 
dp 

dX  (!,».') 

.§£1  Radiation  Boundary  Condition  for  the 
kip)  -  fiy  +  |  yth  propagating  mode 

Figure  6.  Fluid  Domain 


-p(-W) 
=  dp 

dx  (-hyp) 
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wave  equation  numerically  we  must  derive  an  equivalent  finite  difference 
formula,  from  which  we  can  solve  for  the  pressure  for  ail  subsequent  time 
levels  as  follows: 


7.2  a2p  _  a2p  z2p 
' ~dx 2  ' 


f  a/2 


5y2 


(63) 


is  the  two-dimensional  wave  equation  as  derived  in  Chapter  II  Section  A. 
Using  central  difference  approximation  for  all  second  order  partial 
derivatives  the  equivalent  finite  difference  equation  for  Equation  63  is 


kL 

At2 


„n+ 1  ,  „n-ll_  1  [„ri  ,  „n  1,1  f^ri  n  n  ,  „n 

Pi,j  -  2P,,j  +  Pi,j  |  -  pi PM.i  -  2Pi,j  +  Pi-1, j J  +  p[P/,;rt  -  2Pi,i  +  Pi,,- 1 


(64) 


where  h  =  Ax  =  Ay  is  the  step  size  and  At  is  the  increment  in  time.  The 

truncation  error  for  Equation  64  is  0(/i2)  in  space  and  0{At2)  in  time.  Solving 
for  p";+1  explicitly  we  have 


2  Ar 

Letting  p  =  Equation  65  can  be  written  as 


p"}' =2(1~2 p2Yli +p2[^w +phi+p",h  i +p."/-i]-p"7 '• 


The  Von  Neumann  stability  criterion 


P 


(66) 


(67) 


22 


must  be  satisfied  to  ensure  the  stability  of  Equation  66.1 

Special  attention  must  be  paid  when  applying  the  wave  equation  along 
the  boundaries  at  x  -  1  and  -1  for  it  is  here  that  we  make  use  of  the  periodic 
boundary  conditions.  Applying  Equation  66  at  x  =  -1,  i  =  0  and  y  =  jh  (see 
Figure  7)  yields 


This  requires  the  value  of  p  at  the  point  (-1,  jh,  nAt)  which  lies  outside  the 
domain.  By  using  the  periodic  boundary  conditions 

pi  1/  y,  t )  =  pi-1,  y,  t)  (69) 


dp  __  dp 

dx  (+hy,f)  dx 


(70) 


we  can  substitute  the  value  of  (( N-l)h ,  jh,  nAt )  (where  N  is  the  total  number 
of  subdivisions  in  the  x  direction)  for  the  value  at  (~1,jh,  nAt),  allowing  us  to 
evaluate  the  wave  equation  at  the  boundary. 


B.  APPLICATION  OF  THE  RADIATION  BOUNDARY  CONDITION  IN  THE 
NUMERICAL  SCHEME 

As  was  mentioned  at  the  end  of  Chapter  I  (Section  C),  we  apply  a  non¬ 
local  radiation  boundary  condition  (referred  to  as  nlrb )  to  the  fluid  to 
simulate  an  infinite  domain  in  the  positive  y  direction.  Our  domain  is 
truncated  at  jc=1  and  -1,  forcing  our  fluid  domain  to  act  as  a  waveguide.  The 
scattered  pressure  can  be  represented  as  a  series  of  plane  waves  which  take  the 
form 


1For  treatment  of  the  von  Neumann  Stability  Criterion  see  Appendix  A. 
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j[Yk*+Pky-t) 


where 


for  propagating  modes 
for  evanescent  modes 


(71) 


and  yk  =  kn.  We  note  that  when  /3jt  =  - kj  Equation  71  yields 

e-Pky+Knx-t) 

which  is  an  exponentially  decaying  quantity  for  positive  values  of  y.  This  fact 
allows  us  to  neglect  evanescent  modes  when  applying  the  radiation  boundary 
condition,  at  y  =  2. 

The  total  scattered  pressure  is  given  by 

p-  <7» 

k=-°o 


This  is  composed  of  propagating  and  evanescent  modes.  Far  from  the  fluid 
solid  interface  where  only  the  propagating  modes  are  assumed  present  (for 
reasons  given  above)  the  scattered  pressure  is  written 

V=  (73) 

k=-M 


where  M  is  the  total  number  of  modes  (positive  or  negative)  under 
consideration.  At  the  boundary  y  =  2,  we  apply  the  nlrb  operator  (Scandrett 
and  Kriegsmann,  1992,  unpublished  paper). 

=  <74> 


xu  the  individual  modes  of  the  scattered  pressure  of  Equation  73  which  yields 
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BM-  X  |/V(r‘x+fcy'')V  X  ,)ll-  (75) 

fc=-Md^V  J  k=-M  ^dtK  }J 


Equation  75  reduces  to 


M 


b(p)=  YMh-ihW^^- 

k=-M 


(76) 


The  right-hand  side  of  Equation  76  is  identically  zero.  The  boundary  operator 
has  annihilated  the  propagating  modes  and  since  the  evanescent  modes  are 
assumed  to  have  negligible  magnitude  there,  any  scattered  pressure  waves 
reaching  the  boundary  experience  no  reflection,  simulating  an  infinite 
domain  in  the  positive  y  direction. 

To  apply  the  nlrb  operator  at  the  boundary  y  =  2,  we  rewrite  Equation  75  as 


dp 

3y  y=]h 
t=T 


=  0  (77) 


where  we  evaluate  the  pressure  at  a  constant  time  T  and  along  the  boundary 
y  =  Jh  (i.e.  at  y  =  2).  We  incorporate  a k  and  e~iT  as  <ike~iT  and  define  this  to  be  a 


new  constant  ak(T).  Note  that  |Ja^(T)||  =  |fl^|  since  e  ^jj  =  1  for  all  values  of  T. 

<Xk(T)  is  unknown  so  we  must  derive  an  alternative  expression  to  be  able  to 
evaluate  Equation  77.  The  kih  propagating  mode  can  be  written  as 


(78) 


and  when  evaluated  at  y=//i  and  at  f  =  T  and  employing  our  new  constant 
<ty(T)  we  obtain 
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(79) 


t=T 

To  isolate  a]c(T)el^k^h  we  multiply  Equation  79  by  e~lYk*  (we  apply 
orthogonality)  and  integrate  over  the  x  domain  to  obtain 

l 

jp(x,Jh,T)e~tYkXdx  =  2ak(T)el^h  (80) 

-1 

or 

i  1 

ak(T)e^kJh  =  -  J  p($,jh,T)e~lYkSdt,  (81) 

where  £  is  a  dummy  variable  of  integration.  We  substitute  this  value  of 
ak{T)etfikJh  back  into  Equation  77  to  obtain 

+  <82) 

-7  t~T  k=-M  -i 

which  allows  us  to  apply  the  nlrb  at  the  boundary  y  =  2  and  from  Equation  81 
we  will  be  able  to  evaluate  the  amplitudes  of  the  propagating  modes  of  the 
scattered  pressure. 

We  now  proceed  with  deriving  the  finite  difference  approximation  for 
the  radiation  boundary  condition.  Using  a  central  difference  approximation 

for  —  _  «  ,  Equation  82  is  written  as 
dy  y~>n 
3  t=T 
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(83) 


The  trapezoidal  rule  for  integration  is 

fj(x)dx  =  |(/i  +  2/2  +  2/3...+2/„_.,  +  /„)  (84) 


and  is  used  for  the  integral  in  Equation  83,  however  it  can  be  more  compactly 
expressed  as 


where 


1 


r  =  l  or  I 
elsewhere. 


(85) 


where  I  is  the  total  number  of  nodes  in  the  x  direction.  When  substituted 

vu 

into  Equation  83,  using  a  central  difference  approximation  for  (q,  )h,  t). 
Equation  83  can  be  written 


„n  _„n  h*  i 

(86) 

Zk=-M  7=1 


2h 


2  At 


Multiplying  by  ~ jr-,  we  have 


I  M 


m)+2  (87) 


r=l  k=-M 


Define  A  to  be  a  matrix  whose  i,r  entry  is  given  by 

*(.»= 

k=-hi 


(88) 


Upon  substitution  into  Equation  87  we  obtain 
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(89) 


2  At 


=  0. 


We  note  that 


n  n  n  +  1  n- 1 

PiJ+l'PiJ-1  'PrJ  and  PrJ 


are  all  vectors  due  to  the  ] 


subscript  as  described  at  the  beginning  of  this  chapter.  The  radiation 
boundary  condition  and  the  wave  equation  must  both  be  satisfied  at  the 
artificial  boundary  y  =  2.  Applying  either  of  the  two  conditions  will  require 
the  use  of  pseudonodes  which  lie  outside  the  domain.  Through  the 
combination  of  the  two  equations  we  will  be  able  to  eliminate  this 
requirement.  Reproducing  the  two-dimensional  wave  equation  and  the 
radiation  boundary  condition,  (substitute  k  for  all  x  indices  in  the  wave 
equation  and  radiation  boundary  condition,  since  these  are  dummy  indices), 

urg  Hayg 


(91) 


i\0ie  trial  ulna  equations  axt?  u6iu£  6Vdiuat6u  at  y  —  z.  anu  mai  me  vectors 
Pk-\ )  Pk+\  ]  have  a  circular  shift  and  are  of  the  same  dimension  as  the 

rest  of  the  vector  in  Equation  91.  Collecting  like  terms  in  Equations  90  and  91 
we  obtain 


2Pk,}+\~^^Pk,J-l  +pk!jl  +Pk,Jl  +P*-a/ "  ipk.j]~2pk.J  "  0 


kjh 


(92) 
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(93) 


2  At 


2At 


TPll+l^Pl,-i  +  APu-APu-° 


At 


Adding  — j  times  Equation  93  to  Equation  92  yields 


2k 


f 


-2pIh 


f  ry\ 

At2 

kjh 2 


.  „n+l 

+p« 


I  H - 7?  A 


V 


2k 


1 


+  PU 


v 


2fc 


/ 


J 


f  AAt1 
k}h2 


-2 


\.n  df2  n  At  n  _ 

Pk,J  ^Pk-i,l 


(94) 


n+1 


Solving  for  p  ^  j  in  Equation  94  we  obtain 


(  \ 

-l- 

„«+l 

p*,/  = 

7  At  . 

1  H - 7tA 

2kj 

u 

2df2  „ 


df2  „ 


kjh 2  Pk>}~'  +  kjh 2 


2- 


4dr 
kjh 2 


fk.l 


At  n 
+  k2h2Pk^, 


(  \ 

1  — ~~2  A 
2k) 


1C-1 


J 

n  n 


(95) 


The  terms  differenced  in  x,  that  is  Pk-\  j  >  Pk] '  Pk+lJ  can  more  compactly 

expressed  as  Tplj  where  T  is  a  tridiagonal  matrix  with  2--^-  on  the  main 

diagonal  and  on  the  sub  and  super  diagonals.  We  must  also  allow  for 

the  periodic  boundary  conditions  when  constructing  T.  To  do  this  we  replace 

At 2 

the  (N,  2)  and  the  (1,  N-l)  elements  of  T  with  -?-a-  where  T  is  an  N  by  N 

kjh 

matrix.  The  general  form  of  T  can  be  seen  from  Equation  16  Appendix  C. 


Equation  95  is  now  written  as 
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(96) 


„n+l 


f 

-l  - 

r  At 

I  +  ~  A 

2k} 

V  i  J 

2At 


2 TfPh-^Tph- 


kfh 


'  At? 

1 - =rA 

2kf  2 


fk,J 


which  satisfies  both  the  wave  equation  and  the  radiation  condition  at  the 
boundary. 


C  FINITE  DIFFERENCE  APPROXIMATION  FOR  THE  ELASTIC  WAVE 

EQUATION 

To  investigate  the  propagation  of  disturbances  in  an  "I-beam" -shaped 
domain  as  shown  in  Figure  5,  it  will  be  necessary  to  apply  the  elastic  wave 
equation  to  points  in  the  interior  of  the  domain.  (The  boundaries  will  be 
dealt  with  in  a  separate  section.)  By  only  considering  displacements  in  the  x 
(lateral)  and  y  (transverse)  directions,  the  problem  becomes  one  of  plane 
strain  in  two  dimensions.  Reproducing  the  scaled  equations  derived  earlier 
for  motion  in  the  lateral  and  transverse  directions 


1  d2u  1  d2u  f  1 


dzv 


d2U 


kl  dx2  1  k}  ay2  4 

1  d2V  .  1  d2V  .  | 

{ *2 

'  1 

jdxdy  dt2 

d2u  _  d2v 

k2  dx2  kl  dy 2  \ 

k2  jdxdy  dt2 ' 

(97) 


mn\ 

\yo; 


we  use  central  difference  formulas  for  the  partial  derivatives  in  Equation  97 
to  get  the  equivalent  finite  difference  equation 

TjJ  ui+\,j  ~  2ui,j  +  Ui  l,j\ +  ~j2tf  [“JI/+1 "  2Uh  +  “‘V1] 


h\ 
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~  + ^V-). 

=  -T-kw1”2«?»  +  M?  71!  (99) 

Ai  1 '*  '*  '*  • 

The  truncation  error  for  Equation  99  is  0(h2)  in  space  and  OiAt2)  in  time. 
n+ 1 

Solving  for  u  ■  explicitly  in  Equation  99  we  obtain 


W+1  At__ 

M  au2 


[“^+<v]+£[“^'+“"h] 


Ah2  [  k 


■'  1  L,n  ,,n  _,n  ,  ,.n 

i'tf  l  i+l,j+l  ~  i~l'j+1  ~  *+W-i 


.4  »2  ^  i  'N 

-»  1  ,  1  ,!»  -I,”-1 

2  ~  ^  “f*  ^  •  j  ti’  ■  • 

h?-  uf  k }) 


Using  the  same  method  for  Equation  98  we  solve  for  to  obtain 


y”+*  =  ~—(Vn 

klh1 1  H 


K>  .i^hiVuaK^ln) 


.  At2  f  1 _ 1_  /  n  _  n  _un  ...n  ^ 

+  4Ji2[jt£  Jt|j\  ,+1'J+1  ,’~1d+1 


(100) 


„  2dfz  1  1  „  „>! 

+  2 - -S-  +  -S-  Vfi-Vii  . 

I  *2  U  tf  II  ■" 


(101) 


To  ensure  the  stability  of  Equations  100  and  101  the  Von  Neumann  stability 
condition  of 
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At  < 


(102) 


must  be  satisfied.2 


D.  APPLICATION  OF  THE  STRESS-FREE  BOUNDARY  CONDITIONS  TO 
THE  SOLID 

The  boundary  conditions  of  the  two-dimensional  domain  are  broken 


f±n 

down  into  two  major  categories,  those  whose  normal  vector  is  J  ,  where 
txx  =  xXy  =  0,  and  those  whose  normal  vector  is  where  Xyy  =  xXy  =  0. 
These  are  in  turn  divided  into  two  classes.  For  n  =  f  q  j  they  are 

al.  The  boundary  whose  unit  normal  is  fn^| ,  that  is  facing  in  the  positive  x- 

w 

direction,  the  surface  El  in  Figure  5  and 

a2.  The  boundary  whose  unit  normal  is  j^1  j  ,  facing  in  the  negative  x 
direction,  the  surface  DH  in  Figure  5. 

Similarly  for  n  ~  ^  p 

bl.  Boundary  whose  unit  normal  is  f-j  j  ,  the  surfaces  AB,  GH  and  IJ  in 
Figure  5  and 

b2.  Boundary  whose  unit  normal  is  ^  j  ,  the  surfaces  CD,  EF  and  KL  in 
Figure  5. 

The  application  of  the  stress-free  boundary  conditions  for  cases  al  and  bl  is 
discussed  below. 


2For  a  brief  treatment  of  the  Von  Neumann  stability  criteria  see  Appendix 
B. 
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Concentrating  on  displacements  in  the  lateral  ( x )  direction,  a  typical  boundary 
point  as  shown  in  Figure  8  must  satisfy  the  boundary  conditions.  Equations 
103  and  104  and  the  governing  equations  of  motion.  Equations  105  and  106. 

If  we  apply  Equations  105  and  106  at  the  node  (N,y)  in  Figure  8  we  will 

rrt/inir£i  uslivac  fnr  1i^  71^.  71?.  .  .  W?.  .  .  _  ^  T1 H  17?.  _  .  - 

which  lie  outside  of  the  domain  and  are  called  pseudonodes.  To  eliminate 
this  reliance  the  technique  as  developed  by  Ilan  and  Lowenthal  (Ilan,  1976,  pp. 
431-453)  is  followed,  and  is  presented  here. 

The  lateral  displacement  ( u )  at  the  node  (N-l ,/)  in  Figure  8  can  be 
expanded  in  a  Taylor  series  as 


h2  ( -)2  \  higher 
WN-1  /  =  UN  i  ~  h\  T- I  +  T  TT  I  +  order  / 

'1  '1  Ux }frl,j,n  2  )b!,j,n  terms 


(107) 


where  denotes  the  value  of  u  at  the  node  (N-l,;)  and  at  time  level  n. 

Bu  B^u 

Alternate  expressions  for  and  ^2  are  given  by  equations  104  and  105,  we 


have  from  Equation  104 
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From  Equation  105, 
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Thus  the  lower  order  terms  of  Equation  107  can  be  written  as 
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Using  centered  differences  for  ^  an<^  the  following  difference 

formula  for  the  cross  derivative  term 


d2v 


~  ~^h2  (VN,j+l  vN,j- 1  ®N-1J+1+1,N-1,h)  (112) 


3x3y  2 h2 

the  finite  difference  approximation  for  Equation  111  is 
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Cancelling  common  factors  Equation  113  reduces  to 
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Solving  for  u  N  ■  explicitly  we  obtain 
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The  truncation  error  for  Equation  115  is  0(/z3)  (Ilan,  1976,  pp.  431-453).  Using 

the  same  procedure  for  the  transverse  displacement  an  explicit  expression 
n+l  . 

v  „  .•  is  given  as 
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2.  Application  of  the  Stress  Free  Boundary  Condition  for  the  Case  of 


In  the  case  of  b2,  the  governing  equations 


_  du  dv 
tx*~dy  +  dx 


=  0 


(117) 


i  a2u  i  d2u 

_  I  _  -  -  i 

f  1 

l] 

a2z; 

_  d2w 

kl  dx 2  kf  dy2 

1*2 

*?, 

ax3y 

~  at2 

i  a2u  i  a2y 

'  1 

l) 

a2u 

_  d2V 

a  'oJ2  +  z,2  3,  ,2  + 

AJ  A£  1 jif 

72' 

^2 

r  v  m 

«/ 

0/2 

VI- 

(118) 
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(120) 


and  a  typical  boundary  point  is  depicted  in  Figure  9.  Expanding  in  the  vertical 
direction  at  the  node  ( i ,  M+l)  and  at  time  level  n,  and  ignoring  higher  order 
terms 
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using  the  substitution 
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from  Equation  117  and 


from  Equation  119,  Equation  120  now  becomes 
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Figure  9.  Boundary  with  Normal 


dv  d^-u  d^-u 

We  use  central  difference  formulas  for  an£*  for  the  cross 

02  v 

derivative  term  the  following  finite  difference  formula  is  used 
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The  analogous  finite  difference  equation  is  now 
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Cancelling  common  factors  yields 
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and  solving  for  uifM  explicitly,  we  obtain 
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Similarly  for  v  ,•  M  we  have 
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E  FINITE  DIFFERENCE  APPROXIMATIONS  FOR  THE  CORNER  NODES 

The  I  beam  shaped  domain  has  four  270°  corners  which  are  identified  in 
Figure  10  as  1,  2,  3  and  4.  The  treatment  of  these  corners  fails  into  two 
categories,  a)  corners  1,3,  and  b)  corners  2,4.  Within  each  category  the  finite 
difference  formula  applied  is  identical,  the  difference  between  them  comes  in 
the  cross  derivative  terms  of  the  elastic  wave  equation.  Each  category  will  be 
discussed  below.  It  is  important  to  note  that  only  the  governing  equations  of 
motion  are  applied.  The  stress  free  boundary  conditions  are  omitted  due  to 
the  complexity  that  arises  in  trying  to  apply  stress  conditions  at  the  corner 
node.  We  assume  as  in  Fuyuki  and  Matsumoto  that  the  consequences  of 
neglecting  these  boundary  conditions  is  minimal.  (Fuyuki,  1980,  pp.  2051- 
2069) 

Category  a  (Corners  1,3):  An  arbitrary  numerical  mesh  with  Ax  =  Ay  -h 
about  corner  1  is  presented  in  Figure  11.  The  governing  equations  of  motion 
are 
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Figure  10.  The  Corner  Nodes 
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Figure  11.  Comer  Node — Category  a 

In  Equation  124  which  describes  lateral  motion,  the  normal  central  difference 
.  ,  ,  ,  ,  d2u  o2u  B2u 

formulas  are  used  for  the  — r ,  ,  ITT ,  terms  and  for  the  cross  derivative 


d2V 


dt2  '  dx2  '  By2 


term  the  difference  formula  of  Fuyuki  and  Matsumoto  (Fuyuki,  1980, 

pp.  2051-2069)  is  applied  at  the  node  (i,j)  which  is 
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X 

Where  D+  is  the  forward  difference  formula 
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(133) 


and  D  is  the  backward  difference  formula 
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The  resultant  difference  formula  for  the  cross  derivative  term  in  v  at  node  i,j 
is 
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which  is  0(/i2).  The  finite  difference  equivalent  of  Equation  124  is 
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Similarly  for  Equation  125  an  explicit  expression  for  y  is 
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For  Category  b),  the  governing  equations  of  motion  are  the  same.  The  finite 
difference  approximation  of  all  terms  with  the  exception  of  the  cross 
derivative  are  done  in  the  same  manner.  For  the  cross  derivative  term  in  the 
case  of  Equation  124,  the  difference  formula  applied  at  the  node  (i,j)  is 
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Thus  at  that  node  in  Figure  12  is  represented  by 
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The  finite  difference  approximation  to  Equation  124  is  now 
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Figure  12.  Corner  Node — Cateogory  b 
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in  the  same  manner  an  explicit  expression  for  i?--  is 
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F.  BOUNDARY  CONDITIONS  AT  THE  FLUID/SOLID  INTERFACE 


Ihe  boundary  conditions  at  the  fluid/ solid  interface  are  modified  by  the 
introduction  of  a  normal  plane  wave  incident  on  the  surface  of  the  solid.  As 
a  result  the  normal  component  of  stress  is  no  longer  zero.  We  must  allow  for 
the  effect  of  the  scattered  pressure  at  the  interface,  which  is  a  result  of  the 
compliance  of  the  I-beam  structure  and  reflections  of  displacement  waves 


from  internal  boundaries.  This  is  done  by  use  of  the  compatibility  condition 
as  derived  earlier.  Our  necessary  equations  are 
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dp  d2v 


(148) 
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where  t  -  py/ps  and  p5  is  the  scattered  pressure  along  the  boundary  y  ~Kh  (at 
the  interface,  see  Figure  13),  and  at  time  level  n. 

Explicit  expression  for  w"^1  and  are  derived  in  the  same  manner  as 
before  and  they  are 
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Figure  13.  Fluid-Solid  Interface 


Looking  at  the  compatibility  condition.  Equation  142  we  use  central  difference 
approximations  to  obtain 
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assuming  that  it  is  applied  at  the  boundary  y  =  Kh  in  Figure  13, 
Solving  for  we  obtain 
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all  quantities  on  the  right  hand  side  are  known.  We  now  have  a  method  of 
calculating  ViK-\  whic*1  is  required  when  applying  Equation  149  at  the 
boundary  y  -  Kh. 

Solving  for  the  pn+^  term  of  Equation  149  we  have 
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The  components  of  the  right-hand  side  with  the  exception  of  lie  in  the 
interior  or  on  the  boundary  of  the  fluid  domain  and  can  be  evaluated, 
however  pfK^  is  essentially  a  pseudonode  for  the  fluid  and  is  determined 

from  Equation  153  above.  By  this  method  we  have  now  generated  the 
scattered  pressure  waves  caused  by  the  vibration  of  the  solid  at  the  fluid/solid 
surface. 

G.  PROGRAMMING  CONSIDERATIONS 

The  program  for  this  thesis  was  written  entirely  in  Mat  lab  4.0  Beta 

version  for  two  reasons, 

•  Ease  of  programming 

•  The  ability  to  generate  quality  graphics. 

Although  originally  intended  as  a  linear  algebra  toolkit  the  above  features 
have  caused  it  to  be  used  more  and  more  as  a  high  level  programming 
language.  In  our  case  Matlab  was  convenient  since  the  fluid  and  solid 
domains  are  square  matrices,  which  are  easily  manipulated  in  Matlab. 
Updating  values  is  done  somewhat  differently  than  in  FORTRAN  or  C  and  is 
discussed  below. 

Equation  66  of  this  section  updates  the  interior  points  of  the  fluid  domain 
and  is  given  by 


an  equivalent  FORTRAN  statement  might  look  like 
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DO  10  I  =  2,  K 
DO  20  J  =  2 ,  K 


P(I,J,N+1)  «  2..0*  (1.0-2. 0*  (RHO**2)  )*P(I,J,N) 

&  +  (RHO**2)  *  (P  <1-1,  J,N) +P(I+1,  J,N)+P  (I,  J-1,N)  +P  (I,  J+l,  N)  ) 

&  -P (I, J, N-l) ) 

20  CONTINUE 
10  CONTINUE 

Here  each  value  is  updated  individually  by  using  a  double  DO  loop.  In 
Matlab  this  is  done  by  "shifting"  a  grid  the  size  of  the  interior  around  the 
appropriate  matrix  and  weighting  terms.  The  equivalent  code  in  MATLAB 
would  be 

PNEW  (2:K,2:K)  =  2*(l-2*RHOA2)*PCURR(2:K/2:K) 

+(RHOA2)*(PCURR(l:K-l,2:K)+PCURR(3:K+l,2:K)fPCURR(2:K/l:K-l) 

+PCURR(2:K,3:K+l))-POLD(2:K/2:K). 

where  PNEW  contains  the  new  values,  PCURR  the  current  values  and  so  on. 
All  updating  is  done  in  this  manner  eliminating  the  requirement  for 
multiple  do  loops. 
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IV.  NUMERICAL  RESULTS 


In  an  effort  to  verify  our  code  we  checked  the  behavior  of  the  fluid  and 
employed  energy  conservation  methods  to  check  for  the  consistency  of  the 
coupled  domain.  For  the  fluid  waveguide  we  want  to  ensure  that  the 
propagating  modes  behave  as  expected,  that  is,  they  should  not  reflect  from 
the  artificial  boundary.  This  was  done  by  placing  a  driving  force  of  the  form 

p  =  V(ft,y+rnI~')  (155) 


r~5  j  (ou 

where  /?„  -  and  yn  =  nn  and  kf  =  ~,  (Ref.  Equation  24,  Chapter  1, 

Section  A)  at  the  boundary  y  =  0,  which  excited  the  fundamental,  first  and 
second  modes  (r.  -  0,  1,  21-  The  coefficients  An  had  value  1  for  all  n.  The 
amplitude  of  each  mode  was  measured  at  the  boundary.  As  can  be  seen  from 
Figures  14  a,  b,  and  c  the  fundamental  and  first  mode  approach  the  value  1 
with  the  second  mode  showing  the  same  behavior  but  at  a  much  slower  rate. 

To  check  that  the  coupling  of  the  two  domains  was  working  correctly  we 
eliminated  the  cavities  of  the  I-beam  which  reduced  our  problem  to  one 
which  could  be  solved  analytically.  For  a  normally  incident  plane  wave  only 
the  fundamental  mode  was  excited.  Since  the  incident  wave  displaces  the 
solid  only  in  the  y  direction  v  has  no  x  dependence  and  u  is  identically  zero. 
The  values  of  Aq  and  v  are  given  by  (Scandrett,  1992,  interview). 


A)  =C1 


-•( 


i  1  +  e 


-4ikL\] 


and  v  -  e 
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For  kf  =  1  and  ki  =  0.2542,  c\  and  C2  were  given  by  0.0625  +  0.03i  and 
0.0585  -  0.0374/  respectively.  In  this  case  the  magnitude  of  ||Ao|j  of  0.1211 
compared  favorably  with  the  numerical  solution  of  0.1255.  A  major 
discrepancy  lay  in  the  numerical  and  analytical  values  of  v  (transverse 
displacement).  The  average  absolute  value  for  the  numerical  solution  was 
13.12  while  the  average  absolute  value  for  the  analytical  solution  was  0.0965. 
They  exhibited  similar  behavior  but  the  numerical  solution  was  translated  by 
a  constant  term.  We  believe  this  to  be  a  result  of  there  being  no  displacement 
term  to  compensate  for  the  effect  of  imposing  a  periodic  pressure 
instantaneously  in  time  which  has  caused  our  solid  domain  to  drift  or 
displace,  violating  one  of  our  initial  assumptions  (see  Chapter  I,  Section  A). 
To  nullify  this  effect  we  take  the  time  derivative  of  our  steady  state  solution 
for  v  and  then  compare  with  our  analytic  value  as  can  be  seen  in  Figure  15. 
Again  the  numerical  and  analvtical  solutions  give  close  agreement.  A  second 
check  was  to  apply  energy  conservation  methods  to  our  steady  state  solution, 
from  which  it  can  be  shown  (Scandrett,  1992)  that  the  propagating  modes 
must  satisfy 

M  i  } 

X/*«INI  =  -2A)  Re(A))  =  - ^ Im  }Py  =  0Vy  =  0dx '  <156) 

n=-M  -1 

A/f  0 

The  quantities  Mj8n|A„||  and  -2/?0Re(A0)  for  the  various  values  of  kf 

are  listed  in  Table  1.  Considerable  discrepancies  exist  throughout,  which  may 
indicate  either  an  error  in  the  code  or  the  inability  of  our  simulation  to  model 
the  high  frequency  components  which  exist  at  the  interface.  Another  factor 
which  may  contribute  to  the  failure  of  the  integral  is  the  translation  of  v 
which  was  discussed  previously.  With  this  in  mind  we  must  possibly 
consider  the  following  results  as  being  inaccurate. 
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The  first  aspect  we  look  at  is  the  amplitude  of  the  propagating  modes  of 
the  scattered  pressure.  For  kf  =  7,  (this  we  call  our  original  system)  the  time 
history  of  the  amplitudes  for  the  fundamental,  first  and  second  modes  are 
shown  in  Figures  16,  17,  and  18  respectively.  These  indicate  roughly  the  rate 
of  convergence  of  the  numerical  method. 


TABLE  1.  ENERGY  CONSIDERATIONS 


1 

*/ 

XIaAKII2 

-2/%Re(A>) 

0.5065 

-0.1735 

0.022 

0.0003 

1 

0.475 

0.0279 

-0.2138 

1.2516 

0.12965 

0.0006 

0.055 

2.026 

40.5886 

2.135 

3.782 

3.3446 

-0.6093 

0.5380 

-1.3844 

4 

2.6288 

0.2179 

-0.6506 

4.5585 

7.6195 

0.2967 

0.5685 

6.5572 

-5.72 

0.0835 

-0.4053 

7 

4.0212 

0.0644 

0.0667 

7.5 

0.18715 

0.0846 

0.6723 

8 

4.13916 

0.5225 

0.2738 

9 

4.8779 

0.0742 

-1.5974 

To  gain  insight  into  the  characteristics  of  the  solid,  we  treated  the  flange  at 
the  fluid-solid  interface  as  a  thin  plate  and  tried  different  values  of  kf 
corresponding  to  the  resonant  frequencies  of  a  thin  plate  with  prescribed 
boundary  conditions.  The  first  case  was  a  plate  with  the  fixed  (referred  to  as 
FF)  boundary  conditions  given  by 

v(-l)  =  u(l)  =  v"{-\)  =  v"(l)  =  0. 
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The  resonant  frequencies  are 


CjhrPn7 
n  4a^6(l  -  u) 7 


(157) 


(Lalanne,  1982,  p.  103)  where  v  is  Poisson's  ratio,  h  is  the  thickness  of  the 
plate,  cj  is  the  transverse  velocity  of  the  solid  and  a  is  a  scaling  constant  for 
distance. 

For  the  second  case,  damped  (referred  to  as  CC)  boundary  conditions  were 
prescribed,  which  are  given  by 

u(-l)  =  u(l)  =  i/(-l)  =  u'(l)  =  0. 


The  resonant  frequencies  of  this  system  are  given  by 

_  cThX2„ 

wn  —  rrrr 


(158) 


4a^b(l  -  v) 

where  X?  =  22.37,  x\  -  61.67  and  xf  =  120.9 .  (Lalanne,  1982,  p.  103) 

The  purpose  of  this  experiment  is  to  determine  whether  the  I-beam 
exhibits  behavior  similar  to  a  clamped  or  fixed  plate,  since  it  would  be 
advantageous  (numerically)  if  we  could  substitute  a  periodically  placed 
boundary  condition  for  the  remaining  portion  (lower  flange  and  center  spar) 
of  the  I-beam  rather  than  having  to  calculate  finite  difference  approximations 
for  the  entire  I-beam. 

We  do  this  by  plotting  the  amplitudes  versus  the  corresponding  values  of 
kf  for  the  fundamental,  first,  and  second  modes.  The  presence  of  any  peaks 
would  indicate  a  resonant  type  behavior.  If  any  of  these  peaks  corresponds  to 
a  particular  value  of  kf  for  a  clamped  or  fixed  plate  we  say  that  for  that  value 
of  kf  (and  hence  co),  the  I-beam  behaves  in  a  manner  similar  to  a  plate  with 
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clamped  or  fixed  boundary  conditions.  This  gives  us  a  possible  range  for  kf 
values  on  which  to  concentrate  when  looking  for  resonant  frequencies. 

The  values  of  kf  are  given  in  Table  2.  The  quantities  in  parentheses  are 
the  modes  that  propagated  for  a  particular  value  of  kf  which  is  determined  by 
the  maximum  integer  value  n  can  take  such  that  fin  is  real  (Ref.  Equation  71). 


TABLE  2.  kf  VALUES 


SYSTEM 

CATEGORY 

ORIGINAL 

r 

c-c 

F-F 

1 

1(0) 

1.2516(0) 

0.5265  (0) 

2 

4(0,  ±1) 

3.3446  (0,  ±1) 

2.026  (0) 

3 

7  (0,  ±1,  ±2.) 

6.5572  (0,  ±1,±2) 

4.558  (0,  ±1) 

We  divide  our  values  of  kf  into  three  categories  (for  identification 
purposes)  In  the  first  category  we  include  the  values  of  kf  corresponding  to 
the  first  resonant  frequency  of  the  clamped  and  fixed  plates  and  the  values  of 
kf  for  our  original  system  which  allows  only  the  fundamental  mode  to 
propagate.  Ixx  the  second  category  are  values  of  kf  that  correspond  to  the 
second  resonant  frequency  of  the  clamped  and  fixed  plates  as  well  as  the  value 
of  kf  which  allow  only  the  fundamental  and  first  modes  to  propagate.  The 
third  category  corresponds  to  the  third  resonance  of  the  clamped  ard  fixed 
plates  and  the  first  three  modes  in  our  original  system.  Also  included  in  this 
final  category  are  the  values  of  kf  -  7.5,  8  and  9  which  will  allow  us  to  study 
the  behavior  of  the  second  mode  at  higher  frequencies  in  greater  detail. 
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In  Figure  19  we  show  a  plot  of  kf  values  versus  the  amplitude  of  the 
fundamental  mode.  Computed  values  are  shown  in  black.  A  cubic  spline  of 
the  points  involved  is  shown  in  red.  We  must  note  however  that  this  spline 
is  only  a  possible  representation  of  the  behavior  of  the  amplitude  for  different 
values  of  kf,  since  time  constraints  did  not  allow  us  to  conduct  a  more 
comprehensive  set  of  simulations  from  which  we  could  obtain  an  accurate 
picture  of  the  behavior. 

It  can  be  seen  that  for  kf  =  0.5065  and  2.026  (the  points  labeled  FF1  and  FF2, 
the  I-beam  is  exhibiting  a  resonant  type  behavior.  These  values  correspond  to 
the  first  and  second  resonant  frequencies  of  a  fixed  plate.  It  can  also  be  seen 
that  for  values  of  kf  greater  than  3.5  there  is  little  variation  in  the  amplitude 
of  the  fundamental  model  since  we  are  now  past  the  first  cutoff  value  (after 
which  three  modes  propagate)  of  kf  =  k.  This  leads  us  to  believe  that  at  these 
higher  frequencies  most  of  the  activity  takes  place  in  the  higher  modes. 
However  the  total  energy  calculation  shifts  from  fundamental  to  first  and 
back  to  fundamental  as  can  be  seen  in  Table  3. 

In  Figure  20  the  only  resonant  behavior  is  exhibited  for  the  value 
kf  =  3.446  (the  point  labelled  CC2).  This  frequency  is  just  past  the  first  cutoff 
value  in)  and  the  amplitude  of  the  first  mode  is  double  that  of  the 
fundamental  (compared  with  the  corresponding  value  in  Figure  19)  which 
confirms  that  at  higher  frequencies  energy  is  being  propagated  in  the  higher 
modes. 

In  Figure  21  we  plot  the  amplitude  of  the  second  mode  against  kf.  The 
only  resonant  type  behavior  which  exists  here  is  for  the  kf  =  8,  but  since  the 
amplitudes  involved  are  so  small  it  would  be  difficult  to  draw  an  accurate 
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conclusion  about  the  existence  of  a  resonant  frequency.  Another  aspect  we 
investigate  is  the  proportion  of  the  total  energy  carried  by  each  propagating 
mode  which  is  given  by  the  formula 


p=-M 


where  En  is  the  energy  of  the  mode  and  An  the  amplitude  of  the  mode 
(Kinsler,  1982,  p.  110).  These  quantities  are  summarized  in  Table  3. 


TABLE  3.  ENERGY  CONSIDERATIONS 


En(%) 

__ 

*/ 

Eo 

Ei 

E2 

£-1 

E-i 

0.5056 

100 

N/A 

N/A 

N/A 

N/A 

1 

100 

N/A 

N/A 

N/A 

N/A 

1.2516 

100 

N/A 

N/A 

N/A 

N/A 

«  /\n  c 

z.vzo 

r\r\ 

1UU 

XT  /  A 

J.N  /  /-V 

XT  /  A 
i.N  / 

■H 

XT  /  A 

IN  / 

3.3446 

13.04 

43.47 

N/A 

43.47 

N/A 

4 

13.36 

40.82 

N/A 

40.82 

N/A 

4.5535 

7-47 

46.26 

N/A 

46.26 

N/A 

6.5572 

82.43 

2.29 

6.5 

2.29 

6.5 

7 

58.23 

8.08 

12.81 

8.08 

12.81 
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We  can  see  that  the  proportion  of  energy  carried  by  the  fundamental  mode 
varies  from  100%  to  12.04%  as  kf  goes  from  2.026  to  3.3446.  This  can  be 
accounted  for  as  a  redistribution  of  energy  that  takes  place  as  the  frequency  of 
the  waves  passing  through  the  fluid  wave  guide  pass  the  first  cutoff  value  of 
kf  =  tc.  Similar  arguments  hold  for  higher  values  of  kf. 

Finally  we  take  a  brief  look  at  the  shear  strain  field  generated  in  the  solid. 
We  do  this  for  two  reasons,  one  to  check  the  behavior  at  the  corner  nodes  and 
two,  to  find  out  where  the  maximum  shear  strain  occurs. 

When  treating  the  corner  points  we  neglected  applying  the  stress  free 
boundary  conditions  there  assuming  this  would  have  no  adverse  effect  on 
the  behavior  of  the  solid.  Were  this  assumption  invalid  we  would  expect  to 
observe  singularities  or  other  irregular  behavior.  Provided  in  Figures  22,  23 
and  24  are  snapshots  of  the  solid  at  different  time  levels  which  display  no 
unusual  behavior  at  the  corner  nodes,  leading  us  to  accept  the  assumption  as 
valid. 

From  an  engineering  standpoint  we  are  interested  in  where  the 
maximum  shear  strain  occurs  for  possible  failure  analysis.  As  can  be  seen 
from  Figure  25  the  maximum  occurs  along  the  transverse  borders  of  the 
I-beam  cavity.  We  must  note  that  the  amplitudes  and  frequencies  of  the 
acoustic  pressure  waves  involved  in  our  study  are  not  comparable  to  those 
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which  would  cause  permanent  deformation  or  failure,  but  which  could  result 
in  fatigue  cracking  if  subjected  to  prolonged  periods  of  stress. 
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V.  CONCLUSIONS 


Several  problems  arose  in  implementing  the  numerical  scheme.  It  was 
found  that  the  non-local  radiation  boundary  condition  did  not  completely 
annihilate  the  initial  wavefront  and  a  small  amount  of  reflection  took  place, 
but  over  time  the  effect  of  this  reflection  became  negligible.  Evanescent 
modes  did  not  decay  sufficiently  and  were  reflected  at  the  boundary  adding  to 
the  overall  noise  of  the  problem. 

When  calculating  the  amplitudes  of  the  propagating  modes  at  steady  state 
we  applied  the  formula 

1  , 

Ak  =  \p\y=2etknX(ix  (161) 

which  should  yield  consistent  results  for  any  y  values  in  the  fluid  domain.  In 
the  neighborhood  of  the  fluid  /solid  interface  this  was  not  the  case  as  can  be 
seen  from  Figure  26.  We  believe  this  is  due  to  the  presence  of  high  frequency 
components  in  this  region  and  having  too  coarse  a  mesh  to  effectively 
evaluate  the  integral  there. 

With  a  stepsize  h  of  1/40  our  truncation  error  was  on  the  order  of  10~3. 
To  increase  the  accuracy  we  can  perform  Romberg  extrapolation  or  decrease 
the  stepsize,  thereby  increasing  the  dimensions  of  the  matrices  involved.  The 
latter  was  not  an  option  due  to  the  construction  of  the  code  and  machine 
limitations,  in  that  a  simulation  which  involved  10  to  15  thousand  time  steps 
took  from  10  to  12  hours  to  perform.  Increasing  the  size  of  the  matrices 


involved,  thereby  increasing  the  required  time,  would  not  make  it  a  suitable 
code  for  experimentation  and  timely  results. 
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Figure  26.  Amplitude  of  Propagating  Modes  in  the  Y-doinain 


Aspects  of  the  problem  which  deserve  further  study  would  be  the  use  of 
obliquely  incident  waves  vice  normal  as  was  used  here.  Variation  of  the 
cavity  size  and  its  effect  on  the  amplitude  of  the  propagating  modes  should 
also  be  considered.  As  a  simplifying  assumption  the  cavities  were  considered 
to  be  void.  It  would  be  realistic  to  expect  that  they  may  contain  fluid 
(seawater)  or  an  acoustic  dampening  material.  Modifying  the  model  to 
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account  for  such  conditions  warrants  further  study.  A  final  aspect  that  could 
be  considered  is  a  resonance  analysis. 

Throughout  all  of  our  finite  difference  approximations  we  were  able  to 
maintain  second  order  accuracy  in  space  and  time.  We  did  not  have  to  resort 
to  the  use  of  pseudonodes  outside  of  our  fluid-solid  domain.  Finally  the 
possibility  of  being  able  to  establish  an  acoustic  signature  for  a  double-hulied 
structure  could  open  up  a  new  avenue  of  submarine  detection  for  which  this 
thesis  could  be  considered  a  starting  point 
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APPENDIX  A.  VON  NEUMANN  STABILITY  ANALYSIS  FOR  THE  2-D 

SCALAR  WAVE  EQUATION 


The  general  form  of  the  scalar  wave  equation  being  used  is 

1  d2u  d2u  d2u 


cj  8 12  8xz  8y 


.2 


(A-l) 


whose  equivalent  difference  equation,  given  Ax--  Ay  =  h  is 

~cjAt2(U^  ~~ 2u*'i  +  <i')  =  ("f+L;  " 2ui,j  +  +  ^(ui,j+ 1  " Zui,j  +  ui,j- 1  )' 

(A-2) 


The  error  function  takes  the  form 

u  = 


(A3) 


This  is  substituted  into  Equation  A-2  and  common  terms  are  cancelled  to  give 

^2  (€-2+  r')  -  £(«**  +  -2+r'*).  (A-4) 

Using  the  following  identity. 

Jar  ,  -iux 

cos  ccv  -  (A- 5) 

2 

letting  fill  =  yfi  and  defining  p  ^  cAt/h,  Equation  A-4  reduces  to 

4  -  2  +  $  -1  =  p2(4  cos  ph  ~  4).  (A-6) 

Multiplying  Equation  A-6  thiougk  by  £  and  collecting  like  terms  gives 
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42  -  2$(  1  - lp2(cosjip  - 1))  + 1  =  0 


(A-7) 


or  equivalently 


(  2  •  2//J/1V 

l-4p  sm  — - 

'  l  9  J 

\  v  *  J) 


+  1  =  0. 


The  roots  of  this  quadratic  are 


,2  .,;„2 


4l/42  =  1“4P  sin 


For  stability  we  require 


4*1 


which  forces 


1  -  4 p  sin 


•  2|Wf 

in  —  - 

V  2  JJ 


4  <  0. 


Solving  for  p2  in  Equation  A-10  gives 


p2<-  1 


2sin2f~'| 


v  z  / 


which  reduces  to  our  final  stability  criterion  of 

P*  1 


2sin  2<Ph 


or 


(A-8) 


4  .  (A-9) 


(A-10) 


(A-ll) 


(A-12) 
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APPENDIX  B.  VON  NEUMANN  STABILITY  ANALYSIS  FOR  THE  ELASTIC 

WAVE  EQUATION 

In  this  section  wo  give  a  brief  outline  of  the  stability  analysis  required  for 
the  elastic  wave  equation,  whose  general  form  is 


'«  ,  .2  32«  ,  /  J  J.\Z2V  _d2u 

^  T  +  CTr“2+lcL“C7  —— =  ~y 

oa  dyz  v  1  dxdy  dr 


(B-l) 


2  d  V  2  OP 

cf  —  +  ct~  + 


dxz 


V 


(cl-ct) 


d2u  d2v 


dxdy  dt2 


The  finite  difference  approximation  for  Equation  B-l  is 


(B-2) 


,i  -  2uh + “."-l.,) (“”m  -  K + uIh) 

+  tf-W-i)  -  “a  KJ1  -  K,  +'£’)  ®-3) 


We  use  the  following  error  functions 

and  t  =  V{VW* 


(B-4) 


which  are  substituted  into  Equation  B-3,  common  terms  are  cancelled  and 
complex  exponentials  are  gathered  into  trigonometric  quantities  to  give 

-4r  2 |c2 sin2 ^  +  Cj- sin2 ~ jli  -  r2 (c£  - cj )(sin ph sin  yh)V  =  -2  +  £-1j  -U.(B-5) 


Following  the  same  procedure  for  Equation  B-2  we  obtain 
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-4r2^Cj sin2 +  cj sin2 jv  - r2 (c2  - Cj j(sin/Jfr sin yh)U  ~  ^  -2  +  £  Tj- V.(B-6) 


Equations  B-5  and  B-6  can  be  written  in  matrix  form  which  is 


^r2fc2sin2^  +  c£sin2-y-j  -r2^c2  -cjjsin/i/jsin  yh 
-r2(c2  - cj  jsinfihsin  yh  -4r2^c2  sin2  ~ 


2  •  2 

+  cf  sm  — 


fin 


=  yi 


an 


vvy 


(B-7). 


where  A  =  E>-  2  +  h  The  eigenvalues  of  matrix  in  Equation  B-7  are 

r2/c2_c2\  l 

2r2^  +cjj(cosy3/i  +  cos^-2)± — - -[Jcos(/J/i-  yh)  +  cos(ph  +  yli)~-2)2  j2 


they  take  on  a  maximum  value  when  (5h  ~yh  -  n,  and  we  obtain 

Aj  =  A  2  =-4r2jc2  +Cj). 

We  are  now  left  with  the  identity 

£-2  +  £~l  =-4  r2(cj:+c£) 


2 

(B-8). 

(B-9) 


(B-10) 


or 


£2  +  +  <#)- 2^  + 1  =  0. 


(B-ll) 


For  stability  we  require 


which  forces 
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^r2{cl  +Cj-j-2j  -4^0 


(B-12) 


At£ 


h 

tJCL  +  Cj 


(B-13) 


APPENDIX  C  FINITE  DIFFERENCE  FORMULAE  FOR  THE  EQUATION  OF 
MOTION  AND  BOUNDARY  CONDITIONS 


This  appendix  contains  the  finite  difference  approximations  of  all  the 
equations  required  for  the  simulation. 
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Stress-free  Boundary  Conditior 
Surfaces  with  n  =  fjl 
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Surfaces  with  n  =  1 1 
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Corners  2  and  4 
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Fluid 

1)  Wave  Equation 
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APPENDIX  D.  COMPUTER  CODE 


%  function  basel  -  base  for  the  'I  -  beam  '  shaped  domain 
%  written  by:  Lt.  Hugh  He  Bride  USNR 
%  date  :  Apr  92 

%  constructs  a  grid  in  the  shape  of  an  I  beam 

function  (b,rows,  pcoll, pcolr, fill, coll, coir}  »  basel(n) 

%  variables 

%  rows  -  identifies  the  rows  to  be  zeroed  to  form  the  I  beam 

%  coll  -  identifies  the  columns  to  be  zeroed  to  the  left  of  the  center  spar 

%  coir  -  identifies  the  columns  to  be  zeroed  to  the  right  of  the  center  epar 

%  pcolr  -  identifies  the  columns  to  be  zeroed  to  the  right  cf  the  center  spar 
«  including  the  boundary  values. 

%  pcoll  -  Jentifies  the  columns  to  be  zeroed  to  the  left  of  the  center  spar 
%  including  the  boundary  values. 

%  bb  -  building  block  of  the  correct  size 

%  cc  -  dimension  of  the  blocks  to  be  zeroed  out  on  either  side  of  the  center 
%  spar 

t  s  -  variable  to  adjust  the  size  of  the  zero  blocks  and  keep  it  symmetric 
%  fill  -  block  of  zeros  of  appropiate  size  to  fill  in  the  spaces  on 

%  either  side  of  center  spar  i.e.  zeroing  out, at  every  time  step. 

%  *  -  no.  of  divisions  in  the  half  length  of  the  domain. 

%  cr  -  variable  used  to  pick  the  columns  to  the  right  of  the  spar 


%  build  a  grid  of  the  correct  size, 
bb  *■  ones (2*n+l)  ; 


%  the  basic  variables  required  to  build  the  X  -  beam  shaped  grid 
m  “  n*l; 
s  -  log2 (n) ; 
cc  «*  (s-l)*(log2  (n) )  ; 

%  determine  the  row  numbers  to  be  zeroed  out. 
rows-  m-(cc-l) :m+(cc-l) ; 


%  columns  to  the  left  of  the  center  spar,  pcoll  includes  the  boundary  points 
%  coll  does  not. 
coll  »  lscc+1; 
pcoll  “  l:cc+2; 


%  columns  to  the  right  of  the  center  spar,  pcolr  includes  the  boundary  points 
%  coir  does  not. 

cr  ■  2*n+l  -  (cc) ; 
coir  -  cr:2*n+l; 
pcolr  «  cr~l:2*n+l; 

%  construct  the  block  of  zeros  used  to  fill  in  the  spaces  on  either  side 
\  of  the  center  spar  at  every  time  step, 
e  -  size (rows);  e  -  e(2); 
d  size(coll);  d  «  d(2); 
fill  -  teros(e,cc+l) ; 
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bb( rows, coll)  *  fill; 
bb( rows, coir)  -  fill; 

b  **  bb; 
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%  function  bndfnx  -  boundary  facing  in  the  negative  x  direction 
%  i.e  unit  normal  *=  (-1,0) 

%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

% 

%  bndfnx  applies  the  traction  free  boundary  and  the  elastic 
%  wave  equation  along  the  boundary  of  the  I-beam  with  unit  normal 
%  (-1,0)  the  technique  used  is  that  as  developed  by  Ilan  and  Lowenthal 

%  and  is  not  discussed  here. 

%  the  columns  containg  the  nodes  on  the  surface  in  question  and 
%  the  required  neighbours  (those  one  column  in  from  the  surface) 

%  are  picked  off  from  the  current  and  old  values  of  u  and  V 
%  the  shifted  and  weighted  with  constants  from  the  vectors 
%  cunx  (Constants  for  U-values  for  the  Negative  X  boundary) 

%  and  cvnx  (Constants  for  V-values  for  the  Negative  X  boundary) 

%  and  inserted  in  the  correct  position  of  the  updated  u  and  v. 


function  [un  ,vn]  *  bndfnx(uc(vc,uo, vo,un, vn,rows(pcoll,c,d) ; 


%  variables 

%  un  :  updated  values  of  u  vn  :  updated  values  of  v 

%  uc  :  current  values  of  u  vc  :  current  values  of  v 

%  uo  :  old  values  of  u  vo  :  old  values  of  v 

%  rows  :  used  to  identify  the  elements  of  the  matrix  which 

%  are  zero  for  all  times  they  also  contain  the  row  location  of 
%  the  nodes  on  the  boundary. 

%  pcoll  carries  out  the  same  function  as  rows  for 
%  the  columns,  and  the  contains  the  column  location  of  the 
%  of  the  boundary  facing  the  negative  x  direction 

%  c  =  cunx; 

%  d  =  cvnx; 


( i 3 , j  3 ]  =  size (rows); 

mrows  =  (rows(i3)-i  rows  rows( j3)+l] ; 

[ i 4 , j  4 ]  =  size(mrows); 

(i5,js)  =  size(pcoll); 

%  cul  cu2  cvl  cv2  co  cov  contain  the  necessary  u  and  v  values 
%  for  our  calculations 

cu2  =  uc (mrows, pcoll ( j5) +1) ; 
cul  -  uc(mrows,pcoll( j5) )  ; 

cv2  -  vc(mrows, pcoll (j5) +1) ; 
cvl  =  vc (mrows, pcoll( j5) ) ; 

co  =  uo(mrows, pcoll (j5) )  ; 
cov  =  vo (mrows,  pcoll ( j5)  )  ; 


%  the  updated  values  are  calculated 

% 

ucl  =  c(l) *cul (2  :  j4-l)  -  co  ( 2 :  j4-l )  +  c.(2)  *cu2  ( 2  :  j 4-1)  .  . . 
+c(3 ) * (cul ( 3 : j4 )  +  cul ( 1 : j 4  —  2 ) )  .  .  . 
+C(4)*(cv2(l:j4~2)  -  cv2(3:j4))... 
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+c(5) *(cvl(3: j 4 )  -  cvl (1: j4-2) ) ; 


vcl  *  d(l) *cvl (2 : j4-l)  -  cov ( 2 : j4 -1)  +  d (2) *cv2 (2 : j4-l) . . 
+d(3) * (cvl (3 : j4)  +  cvl(l: j4-2) ) . . . 

+d ( 4 ) * (cu2 ( 1 : j4 -2 )  -  cu2 (3 : j4 ) ) . .  , 

+d(5) * (cul (3 : j4)  -  cul(l: j4-2) ) ; 


%  and  put  in  their  proper  place  in  un  and  vn 

pb=  size  (a, rows)  ; 

pbc  -  mrows(2:pb(2)-l) ; 


un(pbc,pcoll (i5) )  =  ucl; 
vn(poc,pcoll( j5) )  =  vcl; 
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%  function  bndfny  -  boundary  facing  in  the  negative  y  direction 
%  i.e  unit  normal  *  (0,-1) 

%  writtan  by:  Lt.  Hugh  Me  bride  OSNR 
%  date  :  Apr  92 

% 

%  bndfny  applies  the  traction  free  boundary  and  the  elastic 
*  wave  equation  along  the  boundary  of  the  I-beam  with  unit  normal 
%  (0,-1)  the  technique  used  is  that  as  developed  by  Han  and  iowenthal 

%  and  is  not  discussed  here. 

%  the  rows  containg  the  nodes  on  the  surface  in  question  and 
%  the  required  neighbours  (those  one  row  in  from  the  surface) 

%  are  picked  off  from  the  current  and  old  values  of  u  and  v 
A  the  shifted  and  weighted  with  constants  from  the  vectors 
%  cuny  (Constants  for  U-values  for  the  Negative  Y  boundary) 

%  and  cvny  (Constants  for  V-valuer  lor  the  Negative  Y  boundary) 

%  and  inserted  in  the  correct  position  of  the  updated  u  and  v. 


function  (un  ,vn]  ■  bndfny (ue,vc,uo, vo,un, vn,rows,pcoll,pcolr,e,d) ; 

%  variables 

%  un  :  updated  values  of  u  vn  :  updated  values  of  v 

%  uc  :  current  values  of  u  vc  :  current  values  of  v 

%  uo  :  old  vale  ;s  of  u  vo  :  old  values  of  v 

%  rows  :  used  to  identify  the  elements  of  the  matrix  which 

%  are  zero  for  ail  times  they  also  contain  the  row  location  of 
%  the  nodes  on  the  boundary. 

%  pcoll  and  pclor  are  both  required  as  there  are  twe  regions, one  on 
%  either  side  of  the  center  spar  of  the  I  -beam  whicn  require 
%  our  attention  and  they  contain  the  location  of  the  nodes  in  question 


%  c  =  cuny; 
%  d  =  cvny; 


(il, jl)  =  size(pcoll)  ; 

[  i2 , j 2  )  =  size(rows); 

[sr  sc)  =  size(uc) ; 

%  rl**  ru*  rv*  and  ro*  pick  off  the  rows  on  either  side  of  the 
%  center  spar  of  the  necessary  u  and  v  values. 

rlul  =  uc(rows(j2)+l, pcoll) ; 
rlu2  =  uc (rows (j2 ) +2 , pcoll) ; 
rui  =  uc(rovs(j2)+2,pcolr) ; 
rul  =  uc(rows( j2)+i,pcolr) ; 

rlvl  *  vc(rows(j2)+l, pcoll) ; 
rlv2  =  vc  (rows  ( j?.)  +2  , pcoll)  ; 
rv2  =  vc(rows( j2)+2,pcolr) ; 
rvi  =  vc(rows ( j2) +1 , pcolr) ; 

rol  =  uo (rows ( j  2 ) +1, pcoll) ; 
ro  =  uo(rows(j2)+l,pcolr) ; 
rolv  =  vo (rows ( j2 ) +1 , pcoll) ; 
rov  =  vo(rows ( j2 ) +1 , pcolr) ; 
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%  the  updated  values  are  calculated 


ul  =  c (l) *rlul (2 : j l-l)  -  rol (2 : j 1-1}  +  c (2) *rlu2 (2 : j 1-1) . . . 
+c(3) * (rlul (3 : jl)  +  rlul (1: j 1—2 ) ) . . . 

+c  (4 )  *  (r.lv2  ( 1 :  j  1-2)  -  rlv2(3:jl))  ... 

+c(5) * (rlvl (3 : jl)  -  rlvlfl: j 1-2 ) )  ; 

ur  -  c(l)*rul(2: jl-1)  -  ro(2:jl-l)  +  C{2 ) *ru2 (2 : j 1-1) .  . . 
+c(3)*(rul(3; jl>  +  rul(l:jl-2)) ... 

+c(4)*(rv2(l: jl-2)  -  rv2 (3 : jl) ) . . . 

+c (5) * ( rvl ( 3 : j 1)  -  rvl <1: jl-2) ) ; 


vl  =  <3 ( 1) *rlvl (2 : j 1-1)  -  rolv(2:jl-l)  +  d(2) *rlv2 (2 : jl-1) . . . 
+d(3)  *  (rl''l  (3 :  jl)  +  rlvl  (1:  jl-2) ) . . . 

+d(4) *(rlu2(l: jl-2)  -  rlu2 (3: jl) ) . . . 

+d  (5) * (rlul (3 : jl)  -  rlul(l: jl-2) )  ; 

vr  =  d(l)*rvi(2; jl-1)  -  rov(2:jl-l)  +  d(2) *rv2 (2: jl-1) . . . 

+d (3) * (rvl (3 : j 1)  +  rvl(l: jl-2) ) . . . 

+d  (4) * (ru2 (1 : j 1-2)  -  ru2  ( 3  :  j l) )  . . . 

+d(5)*(rul(3: jl)  -  rul(l: jl-2) ) ; 


%  and  put  in  their  proper  place  in  un  and  vn 

pb=  size(pcoll) ; 

pbl  =  pcoll (2 :pb (2) -1) ; 

pbr  =  pcolr(2 :pb(2) -1) ; 


un(rows(j2)+l,pbl)  =  ul; 
un ( rows ( j2)+l,pbr)  =  ur; 


vn(rows( j2)+i,pbl)  =  vl; 
vn(rows( j2)+i,pbr)  =  vr; 

%  the  same  procedure  is  reppaated  for  the  'bottom'  of  tha  Z-beam 

un(l/2:sc-l)  =  c (1) *uc (1,2: sc-1 )  -  uo( 1 , 2 : sc-l) . . . 

+  c (2 ) *uc ( 2 , 2 : sc-1 ) . . . 

+c ( 3) * (uc ( 1 , 3 : sc)  +  uc( 1 , 1 : sc-2) ) . . . 

+c (4 ) * ( vc (2 , 1 : sc-2 )  -  vc ( 2 , 3 : sc) ) . . . 

+c ( 5 ) * ( vc ( 1 , 3 : sc)  -  vc ( 1 , 1 : sc-2) )  ; 


vn(  1,2:  sc-1)  —  u  \  i)  "  Vl  ^  1  ,  ^  i  bt’-i  j  —  vo  (l,z:  sc-1)  .  .  . 
+  d (2 ) *vc (2 , 2 : sc-1) . .  . 

+d(3) * (vc(l, 3 :sc)  +  vc( 1 , 1 : sc-2) ) . . . 

+d (4) * (uc (2 , 1 : sc-2 )  -  uc ( 2 , 3 : sc) ) . . . 

+d( 5) * (uc ( l , 3 : sc)  -  uc(l, l:sc— 2) ) ; 
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%  function  bndfpx.  -  boundary  facing  in  the  positive  x  direction 
%  i.e  unit  normal  «=  (1,0) 

%  written  by:  Lt,  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

% 

%  bndfnx  applies  the  traction  free  boundary  and  the  elastic 
%  wave  equation  along  the  boundary  of  the  I-beam  with  unit  normal 
\  (0,1)  the  technique  used  is  that  as  developed  by  Han  ana  Lowenthal 
%  and  is  not  discussed  here. 

%  the  columns  containg  the  nodes  on  the  surface  in  question  and 
%  the  required  neighbours  (those  one  column  in  from  the  surface) 

%  are  picked  off  from  the  current  and  old  values  of  u  and  v 

%  the  shifted  and  weighted  with  constants  from  the  vectors 
%  cupx  (Constants  for  b-values  for  the  Positive  X  boundary) 

%  and  cvpx  (Constants  for  V-valuas  for  the  Positive  X  boundary) 

%  and  inserted  in  the  correct  position  of  the  updated  u  and  v. 

function  [un  ,vn]  =  bndfpx(uc,vc,uo,vo,un,vn,rows,pr,olr,c,d) 

%  variables 

%  un  :  updated  values  of  u  vn  :  updated  values  vh  v 

%  uc  :  current  values  of  u  vc  :  current  valuer,  if  v 

%  uo  :  old  values  of  u  vo  :  old  values  of  v 

%  rows  :  used  to  identify  the  elements  of  the  matrix  which 

%  are  zero  for  all  times  they  also  contain  the  row  location  of 
%  the  nodes  on  the  boundary. 

%  pcolr  carries  out  the  same  function  as  rows  for 
%  the  columns,  and  the  contains  the  column  location  of  the 
%  of  the  boundary  facing  the  positive  x  direction 


%  c  =  cupx; 
%  d  =  cvpx; 


( i 3 , j  3 ]  =  size (rows) ; 

mrows  =  (rows(i3)-l  rows  rows(j3)+l); 

[14. j  4 ]  =  size (mrows); 


%  cul  cu2  cvi  cv2  co  cov  contain  the  necessary  u  and  v  values 
%  for  our  calculations 


cu2  =  uc ( mrows, pcoli (1) -1) ; 
cul  *=  uc  (mrows,  pcolr  (1)  )  ; 

cv2  =  vc(mrows, pcolr (1) -1)  ; 

Ci»  1  a  o rs  1  1  \  V  ■ 

V  1  *  w  y  iii4.  wnM  ^  £> w t/  j.r  y  u.  y  j  j 


co  =  uo(mrows, pcolr (1) ) ; 
cov  =  vo (mrows , pcolr (1) ) ; 

%  the  updated  values  are  calculated 


ucr  =  c ( 1) *cul (2 : j 4 - 1 )  -  co(2:j4-l)  +  c (2) *cu2 (2 : j4-l) . . . 
+c ( 3 ) * (cul ( 3 : j  1 )  +  cul ( 1 : j4-2 ) ) . . . 

+c(4) *(cv2 (1 : j 4 -2 )  -  cv2 ( 3 : j 4 ) ) . . . 

+c  (5)  *  (cvl  (3  :  j  4)  -  cv.l  (1 :  j4-2) )  ; 
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vex 


■=  d'l)  *cvl(2:  j4-lj  -  cov’(2:j4-l)  +  d  ( 2 )  *cv?.  (  2  :  j4~l)  .  .  . 
+d ( 3 ) *(CVl(3:j4)  +  cvl (1: j4-2) ) . , . 

+d{4)*(cu2(l: j4-?.)  -  cu2(3:)4)) ... 

•t-d { 5)  *  (cul (3  :  j4)  -  cul(l:  j 4 -2 )  }  ; 


%  and  put  in  their  proper  place  in  ui.  and  vn 

pb=  5 i ze ( mrows ) ; 

pbc  ■=  mrows{2 :pb (2) -1) ; 


un(pbc,pcolr(l) )  =  ucr; 
vn'pbc,pcolr (I) )  =  ver; 
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% 

% 

X 

% 

% 

% 

% 

l 

% 

% 

% 

% 

% 

* 

% 

% 


function  bndfny  -  boundary  facing  in  the  positive  y  direction 

i.e  unit  normal  =  (0,1) 

written  by:  Lt.  Hugh  Me  Bride  USNR 

date  :  Apr  92 

bndfpy  applies  the  traction  free  boundary  and  the  elastic 
wave  equation  along  the  boundary  of  the  I-beam  with  unit  normal 
(0,1)  the  technique  used  is  that  as  developed  by  Ilan  and  Lowenthal 
and  is  not  discussed  here. 

the  rows  containg  the  nodes  on  the  surface  in  question  and 
the  required  neighbours  (those  one  row  in  from  the  surface) 
are  picked  off  from  the  current  and  old  values  of  u  and  v 
the  shifted  and  weighted  with  constants  from  the  vectors 
cupy  (Constants  for  U-values  for  the  Positive  Y  boundary) 
and  cvpy  (Constants  for  V-values  for  the  Positive  Y  boundary) 
and  inserted  in  the  correct  position  of  the  updated  u  and  v. 


function  [un  ,vn]  =  bndfpy(uc,vc,uo,vo,un,vn,rows,pcoll,pcolr,c,d) ; 

X  variables 

%  un  :  updated  values  of  u  vn  :  updated  values  of  v 

%  uc  :  current  values  of  u  vc  :  current  values  of  v 

%  uo  :  old  values  of  u  vo  :  old  values  of  v 

%  rows  :  used  to  identify  the  elements  of  the  matrix  which 

%  are  zero  for  all  times  they  also  contain  the  row  location  of 
%  the  nodes  on  the  boundary. 

t  pcoll  and  pclor  are  both  required  as  there  are  two  regions, one  on 
X  either  side  of  the  center  spar  of  the  I  -beam  which  require 
%  our  attention  and  they  contain  the  location  of  the  nodes  in  question 


%  c  =  cupy; 

%  d  =  cvpy; 


[il,jl]  -  size (pcoll) ; 

[ i2 , j  2 ]  =  size (rows); 

%  rl**  ru*  rv*  and  ro*  pick  off  the  rows  on  either  side  of  the 
%  center  spar  of  the  necessary  u  and  v  values. 


rlui  =  uc (rows (1) -1 , pcoll) ; 
rlu2  =  uc(rows(l)-2, pcoll) ; 
ru2  =  uc (rows (1) -2, pcolr) ; 
rui  =  uc (rows (l) -1, pcolr) ; 


rlvl  *>  vc(rows(l)-i, pcoll)  ; 
rlv2  =  vc (rows ( 1) -2 , pcoll) ; 
rv2  =  vc(rows(l) -2, pcolr)  ; 
rvi  -  vc(rows(i) -1, pcolr) ; 

roi  =  uo (rows ( 1) -l ,  pcoll )  ; 
ro  =  uo (rows (l) -1 , pcolr)  ; 
rolv  =  vo ( rows ( 1) -i , pcoll )  ; 
rov  =  vo (rows ( 1 ) -1 , pcolr) ; 


94 


%  the  updated  values  are  calculated 


ul  «  c(l) *rlul (2: j 1-1 )  -  rol (2 :  j 1-1)  +  c(2) *rlu2 (2: jl-i) . 
+c(3) * (rlul(3 : jl)  +  rlul (l: jl-2) ) . . . 

+c(4) *(rlv2(l: jl-2)  -  rlv2(3: jl)) ... 

+c(5)*(rlvi(3: jl)  -  rlvl (l: jl-2) ) ; 

ur  ■  c(l) *rul ( 2 : j 1-1)  -  ro(2:jl-l)  +  c(2) *ru2 (2: jl-l) . . . 
+c (3) * (rul (3 : j 1)  +  rul(l: jl-2) )  . .  . 

+c(4) * (rv2(l : jl-2)  -  rv2 ( 3 : j 1) ) . .  . 

+c(5) * (rvi ( 3 : j 1)  -  rvl (l: jl-2) ) ; 


vl  «=  d(i)*rlvi(2: jl-l)  -  rolv(2:ji-l)  +  d(2 ) *rlv2 (2 : j l-i) 
+d(3) *(rlvl(3: jl)  +  rlvl (1 : jl-2) ).. . 

+d(4) *(rlu2 (1: jl-2)  -  rlu2(3:jl)) ... 

+d(5)*(rlul(3: jl)-  rlul ( 1 : j 1-2) ) ; 

vr  =  d ( 1) *rvl ( 2 : jl-l)  -  rov(2:jl-l)  +  d (2) *rv2 (2 : jl-l) . . . 
+d (3) * (rvi (3 : jl)  +  rvl(l: jl-2) ) . . . 

+d(4)*(ru2(l: jl-2)  -  ru2 (3 : jl) ) . . . 

+d(5) * (rul (3 : j 1)  -  rul (1 : jl-2)); 


%  and  put  in  their  proper  place  in  un  and  vn 

pb=  size(pcoll); 

pbl  *=  pcoll  (2  :  pb  (2 )  -1)  ; 

pbr  =  pcolr (2 :pb(2) -1) ; 


un { rows ( 1) -1. pbl)  =  ul; 
un(rows(l) -l,pbr)  =  ur; 

vn ( rows ( 1) -1 , pbl)  »  vl; 
vn(rows(l)-l,pbr)  *  vr; 
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%  function  cmodck  -  checks  the  amplitude  of  the  propagating  mode 
%  written  by:  Lt.  Hugh  Me  Bride  USNF 
%  date  :  Apr  92 

% 

t  emodek  is  the  driver  program  for  the  problem. All  the  vaue  of  the  constants 
$  are  defined  here  and  all  quantities  are  scaled  before  they  are  fed  into  the 

%  decks  are  cleared  before  calculation 
record ( 'erase ' ) 
clear 
clg 

axis ( 'auto' ) 

dx  =  input('step  size  «=  '); 

a  =  input ('  length  scaling  factor  «  '); 

omg  =  input ('  time  scaling  factor  (omega)  =  '); 

%  constants 

t  ct:  tranverse  velocity  of  the  solid  (steel) 

%  cl:  longitudinal  velocity  of  the  solid  (steel) 

%  dnss:  density  of  the  solid  (steel) 

%  cf  :  speed  of  sound  in  fluid  (seawater) 
i  dnsf  :  density  of  fluid  (seawater) 

%  epss:  ratio  ot  fluid  to  solid  densities 

ct  =  3200;  cl  =  5900;  cf  =  1500;  dnsf  =  1026;  dnss  =  7700; 
epss  =  dnsf/dnss; 

%  the  determination  of  the  scaled  variables 
%  dxs  scaled  distance 
%  dts  scaled  time 
dxs  =  dx/a; 

kf  =  omg/cf;kt  =  (omg*a)/ct;  kl  =  (cmg*a)/cl; 

%  dtsl  and  dts2  are  the  stability  criteria  for  the  solid  and  fluid 

%  always  choose  the  minimum 

dtsl  =  (kl*dxs) / (sqrt (1  +  (ct*2/clA2) ) ) ; 

dts2  =  . 5* ( (kf*a*dxs) /sqrt (2) )  ; 

(dtsl  dts2] 

dts  =  input('  desired  time  step  =  '); 
nn  *  input ('  number  of  timesteps  =  '  ); 

%  st  -  when  to  stop  building  the  radiation  boundary  condition 
%  for  in  the  construction  of  the  matrix  A  each  loop  thru  the  construction 
%  cylce  allows  another  mode  to  propagate 


i  =  sqrt(-l) ; 

%  x  a  vector  the  length  of  the  domain  used  in  several  places 
%  i.e.  when  integrating  etc. 

%  mm  -  the  mode  being  checked  0-fundamental  etc. 

I  nt  -  no  of  intervals  in  domain 
%  tt  weighting  factor  for  the  trapezoidal  rule 
x  =  -1 : dxs: 1  ; 

m  «=  size(x)  ;  m  -  m(2)  ;  ml  =  zeros (m)  ; 
mm  «  o 

nt  =  ( l/dxs) ; 
tt  =  (m-l)/2; 
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[un, vn,mn, zO]  «  cvrv4 (ml ,kl , kt ,kf , dts,dxs, x, nn, st , epss, tt ,nt ,mm) ; 

%  un  displacment  in  the  x  dirn 
%  vn  displacment  in  the  y  dirn 
%  mn  displacment  of  the  fluid 

%  zO  vector  containing  the  amplitude  of  the  propagating  mode 
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%  function  coupl  -  couples  the  two  medic  at  the  fluid  solid  interface 
%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

% 

%  coupl  applies  the  traction  boundary  conditions  at  the  fluid  solid 
%  interface  where  the  shear  component  is  zero  but  the  normal  component  is 
4  given  by  tau(yy)  =  -p(total)  .which  causes  the  I  beam  to  deform 
%  and  causing  the  propagation  of  scattered  pressure  waves  to  propagate  in 
%  fluid.  To  prevent  cavitation  at  the  interface  we  apply  the  inviscid 
%  form  of  the  Navier-Stokes  equation  (  or  Euler's  equation  )  at  the 
%  the  boundary 


function  [un  .vn.mnpd]  »  coupl (uc. vc , uo, vo , un, vn, c, d.k, dts , dxs, all , bll , epss) ; 
%  variables 

%  un  :  updated  values  of  u  vn  :  updated  values  of  v 
4  uc  :  current  values  of  u  vc  :  current  values  of  v 
%  uo  :  old  values  of  u  vo  :  old  values  of  v 


%  c  =  cupy; 
%  d  *  cvpy; 


i  =  sqrt(-l)  ; 

(sr  sc]  =  size(uc) ; 

%  the  forcing  function 

time  =  exp (-i*k*dts) 'ones ( 1 , sc) ; 


%  the  u  component  requires  no  modification  and  is  treated  in  the 
%  usual  fashion. 

un (sr , 2 : sc-l)  =  c ( 1) *uc (sr , 2 : sc-1)  -  uo(sr , 2 : sc-1) . . . 

+  c (2) *uc(sr-l , 2  :  sc-1) . .  . 

+  c(3 ) * (uc (sr , 3 : sc)  +  uc(sr, l:sc-2) ) . . . 

+  c(4) * (vc(sr-l, 1 : sc-2)  -  vc (sr-1 , 3 : sc) ) . . . 

+  c ( 5 ) * ( vc ( sr , 3 : sc)  -  vc(sr, l:sc-2) ) ; 

%  the  periodic  boundary  condition  for  u 
un(sr,l)  =  c(l) *uc(sr, 1)  -  uo(sr.l)... 

+  c(2) *uc ( sr-1 , 1 ) . . . 

+  c ( 3 ) * (uc (sr ,  2 )  +  uc(sr.sc-l) ) . . . 

+  c (4 ) * ( vc (sr-1 , sc-1 )  -  vc (sr-1 , 2 ) ) . . . 

+  c (5) * (vc (sr , 2)  -  vc(sr.sc-l) ) ; 

un(sr.sc)  =  un (sr, 1) ; 


%  the  normal  component  of  stress  plus  the  incident, ref lected 
%  and  scattered  pressures 

vn (sr , 2 : sc-1)  =  d ( 1 ) *vc (sr , 2 : sc-1)  -  vo (sr , 2 : sc-1) . . . 

+  d  (2 ) *vc (sr-1 , 2 : sc-1) . . . 

+d ( 3) * ( vc (sr ,  3  :  sc)  +  vc  (sr , 1 : sc-2 ) ) . . . 

+d (4 ) * (uc ( sr-1 , 1 : sc-2 )  -  uc (sr-1 , 3 : sc) ) . . . 

+d ( 5) * (uc (sr , 3 : sc)  -  uc (sr , 1 : sc-2 ) ) . . . 

+  ( (2 *dts~2) /dxs) * ( 2*epss*time (2 : sc-1 )  +  epss*al 1 ( 2 : sc-1) ) ; 


%  the  periodic  boundary  condition  for  v 
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vn(sr , 1)  -  d ( 1 ) *vc (sr , 1)  -  vo(sr,l}... 

■*  d(2)  *vc (sr-i ,  1)  .  . . 

+d(3) * (vc(sr, 2)  +  vc(sr,sc--l) )  . . . 

+d (4 ) * (uc (sr-l , se-1)  -  uc(sr-l,2) ) . ,  . 

■+d(5)  *  (uc(sr ,  2)  -  uc ( sr ,  sc-1 ) )  .  .  . 

+  ( (2*dts',2) /dxs)  *  (2*cpss*tirce  (1)  +  epss*all (1, l) ) ; 


vn(sr,sc)  «=  vn(sr,l); 

%  the  compatabi lity  condition. 

anpd  *=  (- (2*dxs) / (dts~2) ) * ( vn(sr, : ) -  2*vc(sr,:)  +vo{sr,:))  +  bii(i, 
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%  function  cwv4  -  coupled  waves  version  4 
%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

t 

\  cwv4  is  the  main  program  which  couples  the  behaviour  of  the  'I-beam'  shaped 
%  solid  medium  with  the  fluid  medium.  The  elastic  wave  equation  and  the  boundary 
%  conditions  , periodic  and  traction  free  are  are  satisfied  for  the  solid,  while 
*  the  scalar  wave  equation  and  the  periodic  and  radiating  boundary  conditions  ar 
t  applied  to  the  fluid. 

% 

%  The  general  steps  of  the  program  are  as  follows 

% 

%  1.  The  basic  parameters  are  determined  ,that  is  the  size  of  the  domains 

%  as  passed  to  it  by  the  driver  program  cmodck.m 

I 

%  2.  All  the  global  variables  are  calculated  for  both  the  fluid  and  the  solid 

%  including  the  matrices  required  for  the  radiation  boundary  condition. 

% 

%  3.  vl  and  v3  are  vectors  used  to  find  the  amplitude,  vl  is  the 

I  weighting  vector  1/2  ill....  l  1/2  for  the  trapezoidal  rule 
%  and  v3  =  exp(i*n*pi*x)  ,  the  othhogonal  vector. 

% 

%  4.  The  I  beam  is  built  by  basel 

% 

%  5.  The  initial  values  of  the  u  and  v  for  the  solid  (un,vn,uc, vc,uo,vo) 

%  and  m  (inn  me  and  mold)  for  the  fluid  are  set  to  zero. 

% 

%  6,  The  elastic  wave  equation  and  the  periodic  boundary  conditions 

%  are  applied  to  u  and  v 

% 

%  7.  The  boundary  conditions  for  the  solid  are  applied 

% 

%  8.  The  fluid  and  solid  are  coupled 


10. 


11. 


% 

% 

% 

% 

% 

% 

% 

%  12 
% 

% 

% 

% 


9  . 


The  freespace  portions  of  the  1  beam  are  zeroed  cut 

The  scalar  wave  equation  and  the  periodic  boundary  conditions  for 
the  fluid  are  applied. 

The  values  of  the  amplitude  are  calculated  and  accumulated 

u,  v  and  m  are  updated  ,that  is  the  new  values  become 

the  current  values  and  the  current  values  become  the  old  values. 


function  [v,un, vn,mn, z2]  =  cwv4 (m,kl,kt,kf ,dts,dxs,x,nn,st,epss,tt,nt,mm) ; 
%  variables 

%  un  :  updated  values  of  u  vn  :  updated  valuer,  of  v 

%  uc  :  current  values  of  u  vc  :  current  values  of  v 

%  uo  :  old  values  of  u  vc  :  old  values  of  v 

%  rows  :  used  to  identify  the  elements  of  the  matrix  which 

%  are  zero  for  all  times  they  also  contain  the  row  location  of 
%  the  nodes  ori  the  boundary. 

%  pcoll  and  pclor  are  both  required  as  there  are  two  regions, one  on 
%  either  side  of  the  center  spar  of  the  I  -beam  which  require 
1  our  attention  and  they  contain  the  location  of  the  nodes  in  question 
%  kl  -  scaled  longitudinal  speed 
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%  kt  -  scaled  transverse  speed 
%  dxs  -  scaled  spacing 
%  dts  -  scaled  time  step 
%  x  -  vector  of  length  of  the  interval 
%  nn  -  no  of  tisiesteps 

%  st  -  stopping  criteria  for  radiation  boundary  condition 
%  epss  -  ratio  of  solid  density  to  fluid  density 
V  nt  -  no  of  intervals  in  domain 

%  mm  -  prpoagating  mode  under  investigation  ,mm  =  0  fundamental 
%  mm  =  1  first  and  so  forth 
%  z?  -  amplitude  of  propagating  mode 


t  Step  1  :  r,c  the  dimensions  of  the  domain,  r-1  ,c-l  and  n-l 
%  the  dimensions  of  the  interior 

[r  c]  =  size (m) ; 
n  «•  r-l;  r  *=  r-l;c  =  c-l; 
axis( 'xy' ) 

%  Step  2:  The  global  variables 


%  For  the  solid 

rho  *  rhov (kl ,kt , dxs,dts) ; 
cunx  «=  ucnx  (kl ,  kt ,  dxs,  dts)  ; 
cuny  =  ucny(kl,'xt,dxs,dts)  ; 
cupx  =  ucpx(kl,kt,dxs,dts) ; 
cupy  *=  ucpy(kl.kt,dxs,dts)  ; 

%  For  the  fluid 

rhof  *=(dts*2)  /  (kf *2*dxs*2)  ; 
rcl  -  (2*dts*2)/(kf*2*dxs'2) 


ccs  =  (2  -  2*rho(l) -2*rho(2) ) ; 
cvnx  -  vcnx (kl , kt , dxs , dts) ; 
cvny  =  vcny(kl,kt,dxs,dts) ; 
cvpx  =  vcpx(kl,kt,dxs,dts) ; 
cvpy  -  vcpy (kl,kt,dxs,dts) ; 


ccf  =  (2  -  4*rhof) ; 


X  The  matrices  for  the  radiation  boundary  condition  allowing  the  propagating 
%  modes  to  pass  through  the  artificial  boundary. 


anew  =  zeros(c+l); 
for  pm  =  0:st 

acurr  =  rbc(r+l,c+l,pm,dxs,kf ) ; 
anew  «  anew+acurr; 

end 

ml  *  (eye(c+l)  +  (dts/ (2*kf *2) ) *anew) ;  ml  =  inv(ml); 
m2  =  (eye(c+l)  -  (dts/ (2*kf *2) ) *anew) ; 
t  =  trid(dxs,dts,c-t-l,kf )  ; 

%  Step  3:  vl  and  v3  calculated  so  as  to  be  able  to  determine  the  amplitude 
vl  =  del2 (2*nt) ; 
v3  =  exp(i*mm*pi*x' ) ; 

%  Step  4:  basel  builds  the  I  -beam  pclor,pcoll  , coll, coir  and  fill 
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%  allow  us  to  identify  the  boundaries  of  the  solid  domain  and  the 
%  corner  nodes. 

(b,rows,  pcoll,pcolr,  fill,  coll,  coir]  =  basel(n/2); 


%  Step  S:  All  values  for  the  fluid  and  the  solid  are  initially 
t  set  to  zero. 

un  -  zeros (size (m) ) ;  uc  =  un;  uo  ■  un; 
vn  *=  zeros(size(m) ) ;  vc  »  vn;  vo  -  vn; 
urn  «  zeros(size(u) ) ;  me  =  mn  ;mold  =  inn; 


for  k  =  linn 

%  Step  6:  The  elastic  wave  equation  is  applied  to  the  interior 
%  of  the  solid 

un(2:n,2:n)  =  ccs*uc(2 in, 2  in)  -  uo (2 : n, 2 i n) . . . 

+  rho(2)  *(uc(im-l,2in)  +  uc(3:n+l,2:n) ) . . . 

+  rho( 1) * (uc(2 in, li n-l)  +  uc(2:n,3:n+l) ) . . . 

+  rho (3) * (vc (1 1 n-l, li n-l)  +  vc(3:n+l,3:n+l) ) . . . 

-  rho(3)  *  ( vc(  1 :  n-l,  3  in+1)  +  vc  ( 3  :n+l,  1  :n-l) )  ; 


vn(2:n,2:n)  “  ccs*vc (2 : n, 2 i n)  -  vo(2 : n, 2 i n) . . . 

+  rho(l)*(vc(i:n-l,2:n)  +  vc(3 in+1 , 2  in) ) . . . 

+  rho (2 ) * (vc (z : r, 1 ; n-l)  +■  vc(2:n,3:n+l) ) . . . 

+  rho( 3) * (uc( 1 :n-l, li n-l)  +  uc(3 in+1, 3 in+1) ) . . . 
-  rho ( 3) * (uc ( 1 i n-l, 3 :n+l)  +  uc(3in+l, lin-I) )  ; 


%  The  periodic  boundary  conditions  for  the  solid 

un(2:n,l)  *  ccs*uc  (2  :  n,  1)  -  uo(2;n,l)... 

+  rho(2)  *  (uc(  i  m-1, 1 )  +  uc(3m+l,l) )  . . . 

■*  rho  ( 1)  *  (uc  (2  i  n,  2 )  +  uc (2  i  n,  n)  )  .  .  . 

+  rho ( 3 ) * (vc(l :n-l , n)  +  vc ( 3 ; n+1 , 2) ) .  .  . 

-  rho(3) * (vc(n n-l, 2)  +  vc(3in+l,n))  ; 

un(2:n,n+l)  =  un(2:n,l); 


vn(2m,l)  =  ccs*vc(2.-n,  l)  -  vo(2:n,l)... 

+  rho(i) * (vc(iin-i, l)  +  vc(3 in+i, 1) ) . . . 

+  rho(2)  *  (vc(2  in,  2)  +  vc(2m,n) ) . . . 

+  rho ( 3) * (uc ( 1 i n-l, n)  +  uc(3:n+l,2) ) . . . 
-  rho (3) * (uc (l sn-l,2)  +  uc(3in+l,n))  ; 

vn(2:n,n+l)  »  vn(2:n,l); 


t  Step  7 i  The  boundaries  of  the  I  beam  and  including  the  corner  nodes 
%  are  treated. 
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(un  ,vn]  =  bndfnx(uc,vc,uo, vo,un,vn, rows , pcoll , cunx, cvnx) ; 

[un  ,vn]  =  bndfny(uc, vc.uo, vo, un,vn,rows, pcoll, pcolr, cuny, cvny) ; 
[un  ,vn]  ■=  lcorny(uc,vc,uo,vo,un,vn,rows,pcoll,cuny,cvny) ; 

[un  ,vn)  =  bndr'pxfuc,  vc,uo, vo, un,vn, rows, pcolr, cupx.cvpx)  ; 

[un  ,vn]  *  bndfpy(uc,vc,uo, vo,un,vn,rows, pcoll, pcolr. cupy,cvpy) ; 
[un  ,vn]  =  Icorpy (uc, vc,uo, vo,un,vn, rows, pcoll, cupy,cvpy) ; 
[un,vn]  =  tcorl3 (uc, vc, uo,vo,un,vn, rows, pcoll,  pcolr, rho) ; 

[uri,vn]  -  tcor24(uc, vc,uo,vo,un,vn, rows, pcoll, pcolr, rho) ; 


%  Step  8  :  The  solid  and  the  fluid  are  coupled  (Note  the  is  where  the  forcing 
%  function  of  the  problem  is  contained  ) 

all  =  xbc(1  bll  *=  me ( 2 , : } ; 

[un  , vn,mnpd]  =  coupl (uc, vc, uo, vo, un, vn, cupy ,cvpy , k,dts,dxs, all , bll,epss) ; 


%  Spep  9:  Both  sides  of  the  center  spar  for  all  values  of  u  and  v  are  zeroed 
%  out  so  as  to  prevent  pollution 


un(rows,coll) 
uc(rows,coll) 
un(rcws,c.olr) 
uc(rows,colr) 
uo(rows,coll ) 
uo( rows, coir) 


=  fill;  vn{rows, coll) 
=  fill;  vc(rows,coll) 
*  fill;  vn(rows,colr) 
=  fill;  vc (rows, coir) 
«  fill;  vo (rows, coll) 
-  fill;  vo (rows, coir) 


fill; 

fill; 

fill; 

fill; 

fill; 

fill; 


%  Step  10:  The  scalar  wave  equation  for  the  fluid 


%  First  the  fluid/solid  interface 

mn(l,2:c)  =  rhof* (mnpd(2 :c)  +  mc(2,2:c) . . . 
+mc(l, l:c~l)  +  mc(l,3:c+l))  +ccf *mc( 1 , 2 : c) . . . 
-mold(l, 2 :c) ; 

%  and  it's  periodic  boundary  condition 

mn(l,l)  --  rhof *(mnpd(  1,1)  +  mc(2,l)... 
+mc(l,c)  +  mc(l,3))  +ccf *mc ( 1 , 1) . . . 

-moid(l,  1)  ; 

mn(l,c+l)  *  inn  (1,1); 


%  The  interior  points  of  the  fluid 

mn(2:r,2:c)  =  rhof  *  (mc(  1  :r-l ,  2  :  c)  -t  me  (3  :  r+1 , 2  :  -.)  .  .  . 
+mc(2  :  r,  .1  :c-l)  +  me (2 :i , 3 : c+l) )  +ccf  *mc  ( 2  :  r ,  2  :  c)  .  .  . 
-mold(2:r,2:c) ; 


%  The  radiation  boundary  condition 

%  Note:  We  are  required  to  multiply  a  matrix  by  a  row  vector, but  to  do  this 
%  it  must  be  transposed  to  a  column  .  Matlab  takes  the  Hermitian  transpose 
1  by  default  ,  so  to  ensure  the  correct  signs  we  must  void  this  effect  by 
%  taking  the  conjugate  of  the  transpose  before  we  do  our  calculations. 

%  The  process  is  then  reversed  so  as  the  updated  value  has  the  correct  dimension 
% 
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si  «  conj (nc(r, : )  ' )  ; 
s2  «  conj (mc(r+l, : ) ' ) ; 

S3  -=  conj  (mold(r+i,  : )  ' )  ; 
int  «=  ml*(rcl*sl  +  t*s2  -  m2*s3); 
mn(r+l,:)  *  conj (int'); 


%  the  periodic  boundary  condition  for  the  artificial  boundary 

mn(2:r,l)  =  rhof* (mc(l :r-l, 1)  +  mc(2:r,2)  +  ... 
jBC(3:r+l,l)+  mc(2:r,c))  +  ccf*mc(2:r, 1)  -  raold(2:r,l); 

mn(2:r,c+l)  »  mn(2:r,l); 


%  Step  11:  Calulates  amd  accumulates  the  values  cf  the  amplitude 
aa  “  abs ( (mn (r, : ) . *vl) *v3) ; 
z2  =  [ z2  aa] ; 

%  Step  12;  The  values  for  u,v  and  m  are  updated. 

uo  -  uc;  uc  =  un; 

vo  =  vc;  vc  »  vn; 

mold  *=  me;  me  »  mn; 


end 
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%  function  del2  -  del  funtion  used  for  thr  trapezoidal  rule 
%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

\  bulids  a  vector  of  appropiate  length  ured  to  weight  the  elements  of  the 
%  quantity  being  integrated 

function  d  «  del2(n) 
a  -  ones(l,n+l) ; 
a(l)  -  .5;  a (n+1)  «=  .5; 
a2  «  ones(l,n+l) /n; 

d  =  a.*a2; 
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t  function  lcorny  -  left  corner  facing  in  the  negative  y  direction 
%  i.e  unit  normal  «=  (o,-i) 

%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

% 

% 

% 

%  lcorny  perforins  the  same  function  as  bndfny  except  it  only  operates  on 
%  4  points  ,  those  which  lie  at  the  extremities  of  the  domain,  the  points 
%  as  if  fig(  )  Periodicity  is  used  for  pseudonodes  which  lie  outside  the  domain. 
4  This  requires  us  to  pick  out  each  nodes  individually  (9  for  each  4  points  , 

%  making  a  total  of  36  )  they  are  weighted  in  the  same  fashion 
%  as  in  bndfny  and  the  updated  values  for  u  and  v  are  calculated  for  both  sides 

function  [un  ,vn]  =  lcorny (uc,vc,uo, vo,un,vn,rov’s,pcoll,c,d) ; 

%  un  :  updated  values  of  u  vn  :  updated  values  of  v 

%  uc  :  current  values  of  u  vc  :  current  values  of  v 

%  uo  :  old  values  of  u  vo  :  old  values  of  v 

%  rows  :  used  to  identify  the  elements  of  the  matrix  which 

%  are  zero  for  all  times  they  also  contain  the  row  location  of 
%  the  nodes  on  the  boundary. 

%  pcoll  allows  us  to  pick  out  the  column  location  of  those  nodes 
\  to  the  left  of  centerline  and  using  the  periodicity  we  substitute 
%  this  value  into  the  corresponding  point  on  the  right  of 
%  centerline. 


%  c  =  cuny; 
t  d  *=  cvny; 


fil.jl]  ■  size(pcoll) ; 
[i2, j  2  ]  =  size(row5) ; 
[i3, j 3  J  =  size (uc) ; 


4  the  updated  values  of  u  to  the  left  of  centerline  are 
%  calculated 

ucc  =  uc(rows( j2) +1, pcoll (1) ) ; 
uccl  =  uc (rows ( j2) +2 , pcoll ( 1 )) ; 
ucr  =  uc(rows(j2)+l,pcoll(2) ) ; 
ucl  -  uc(rows( j2)+l,  j3-l)  ; 

ver  =  vc(rows( j 2 ) +1 , pcoll  ( 2)  )  ; 
vru  -  vc (rows (j2 )+?, pcoll (2) )  ; 
vlu  =  vc(rows( }2)+2,  j3-l)  ; 
vcl  *=  vc(rows(  j2)+l,  j3-l)  ; 

uol  =  uo (rows ( j2 ) +1 , pcoll  ( 1)  )  ; 

ul  =  c(l)*ucc  -  uol  +  c(2)*uccl  +c(3)*(ucr  +  ucl)... 
+c(4)*(vlu  -  vru)  +  c(5)*(vcr-  vcl); 


uccl  =  uc (1 , pcoll (1) ) ; 
ucul  ~  uc(2 , pcoll ( 1)  )  ; 


ucrl  -  uc(l,pcoll (2) ) ; 
ucll  -  uc(l, j 3 — 1 )  ; 

vcrl  =  vc(l,pcoll (2) )  ; 
vrul  “  vc ( 2 , pcoll ( 2 ) )  ; 
vlul  ■  vc(2< j3-l)  ; 
veil  *=  vc(l, j3-l)  ; 

uoll  «  uo  ( 1 , pcoll ( 1) ) / 


ull  «  c(l)*uccl  -  uoll  +  c(2)*ucul  +c(3)*(ucrl  +  ucll)... 
+c(4)*(vlul  -  vrul)  +  c(5)*(vcrl-  veil); 


%  the  updated  values  of  v  to  the  left  of  centerline  are 
%  calculated 


vucc  «=  vc  (rows  (j2)  +1 ,  pcoll  (1) )  ; 
vuccl  <=  vc(rows(j2)+2,pcoll(l) )  ; 
vucr  «  vc(rows( j2) +1, pcoll (2) )  ; 
vucl  =  vc (rows ( j 2 ) +1, j3-l)  ; 

wer  *  uc  (rows  ( j2)  +1 , pcoll  (2) )  ; 
vvru  =  uc(rows(j2)+2,pcoll(2) ) ; 
wlu  =  uc  (rows  ( j 2 )  +2  ,  j3-l)  ; 
wcl  =  uc(rows(j2)+l,  j  3—1)  ; 

vuol  =  vo(rows(j2)+l,pcoll(l) ) ; 

vul  =  d(l)*vucc  -  vuol  +  d(2)*vuccl  +  d(3)*(vucr  +  vucl)... 
+d(4)*(wlu  -  wru)  +  d(5)*(vvcr-  wcl); 


vuccl  =  vc (l, pcoll ( 1) )  ; 
vucul  =  vc ( 2 , pcol 1(1)) ; 
vucrl  =  vc(l,pccll (2) ) ; 
vucll  =  vc(l, j 3  —  l )  ; 

wcrl  =  uc(l, pcoll  (2) )  ; 
vvrul  =  uc (2, pcoll (2) )  ; 
wlul  =  uc  (2 ,  j3-l)  ; 
well  •  uc(i,  j3-l)  ; 

vuoll  =  vo(l,pcoll(l) )  ; 


vull  =  d(l)*vucci  -  vuoll  +  d(2)*vucul  +d(3)* (vucrl  +  vucll)... 
+d(4)*(vvlul  -  vvrul)  +  d(5)*(wcrl~  well)  ; 

%  the  values  put  in  the  correct  positions 

%  and  since  we  have  periodic  boundary  conditions,  we  insert  the  left 
%  hand  value  into  the  cort espondding  right  hand  position 

un(rows(j2)+l, pcoll  (1)  ,  *»  ul; 
un(rows ( j  2 ) +1, j3)  =  ul; 


% 
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vn 'rows (j2)+l,pcol 1(1) )  «=  vul; 
vn(rows(l2)+l, j 3 )  -  vul; 


un(l,pcoll (l) )  -  ull; 
un ( 1 , j  3 )  -  ull; 

vn(l,pcoll(l) )  -  vull; 
vn(l,j3)  *  vull; 
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%  function  lcorpy  -  left  corner  facing  in  the  positive  y  direction 
%  i.e  unit  normal  -  (o,i) 

%  written  by:  Lt.  Hugh  Me  Bride  USNR 
\  date  :  Apr  92 

% 

% 

\ 

%  lcorpy  performs  the  same  function  as  br.dfpy  except  it  only  operates  en 
%  2  points  ,  those  which  lie  at  the  extremities  of  the  domain,  the  points 
%  as  in  fig(  )  Periodicity  is  used  for  pseudonodes  which  lie  outside  the  domain. 
%  This  requires  us  to  pick  out  each  nodes  individually  (9  for  each  2  points  , 

%  making  a  total  of  18  )  they  are  weighted  in  the  same  fashion 
%  as  in  bndfpy  and  the  updated  values  for  u  and  v  are  calculated  for  both  sides 


function  [un  /Vn]  =  lcorpy (uc,vc,uo, vo, un,vn, rows, pcoll,  c, d) ; 

%  un  :  updated  values  of  u  vri  :  updated  values  of  v 

%  uc  :  current  values  of  u  vc  :  current  values  of  v 

%  uo  :  old  values  of  u  vo  :  old  values  of  v 

%  rows  :  used  to  identify  the  elements  of  the  matrix  which 

*  are  zero  for  all  times  they  also  contain  the  row  location  of 
%  the  nodes  on  the  boundary. 

%  pcoll  allows  us  to  pick  out  the  column  location  of  those  nodes 
%  to  the  left  of  centerline  and  using  the  periodicity  we  substitute 
%  this  value  into  the  corresponding  point  on  the  right  of 
%  centerline. 


%  c  =  cupy; 

%  d  =  cvpy; 


%  there  are  now  only  two  values  wnicn  need  to  be  calculated 
%  as  the  other  boundary  facing  in  the  positive 
%  y  direction  is  at  the  fluid/solid  interface  and  requires 
%  a  special  treatment 

[ii,jl]  =  size(pcoll); 

[  j.  2  ,  j  2  ]  =  size(rows); 

[ i3 , j  3 ]  =  size(uc) ; 


%  the  updated  u  value 

ucc  =  uc ( rows ( 1 ) - 1 , pco 1 1 ( 1 ) ) ; 
ucu  =  uc(rows(l)-2,pcoll(l) ) ; 
ucr  =  uc (rows (1) -1 , pcoll (2) ) ; 
uc 1  =  uc ( rows (1)-I,j3-1) ; 

ver  =  vc(rows ,1) -l, pcoll (2) ) ; 
vru  =  vc(rows ( 1) -2 , pcoll (2) > ; 
vlu  =  vc(rows(i) -2 , j3-i) ; 
vcl  =  vc(rows(l)-l, j  3  —  1 ) ; 

uol  =  uo (rows (1) -1, pcoll (1) ) ; 

ul  =  c(l)*ucc  -  uol  +  c(2)*ucu  +c(3)*(ucr  +  ucl)... 
4c(4)*(vlu  -  vru)  +  c(5)*(vcr-  vcl); 


%  the  updated  v  value 


vucc  «  vc(rows(l)-l,pcoll(l) ) ; 
vucu  -  vc(rovs(l)-2,pcoll(i) )  ; 
vucr  ^  vc(rowe (1) -l,pcoll (2) ) ; 
vucl  -  vc(rovs(l)-l,  j3-l)  ; 

wcr  *  uc(rows(l)-i,pcoIl(2) ) ; 
wru  «  uc(rows(l)-2,pcoll(2) ) ; 
vvlu  «  uc(rows (l) -2, j3-l) ; 
wcl  »  uc(rows(l) -1, j3-l)  ; 

vuol  *  vofrows (l) -i,pcoll (1) ) ; 

vul  «  d(l)*vucc  -  vuol  +  d(2)*vucu  +  d(3)*(vucr  +  vucl)... 
+d(4)*(wlu  -  wru)  +  d(5)*(wcr-  wcl); 


%  using  periodicity  we  update  the  values  to  the  left  and  right  of 
%  the  center  spar 

un(rows(l)-l,pcoll(l) )  =  ul; 
un(rows(l)-l,  j3)  ■=  ul; 
vn(rows(l)-l,pcoll(l) )  *=  vul; 
vn(rows(i)-l, j3)  =  vul; 
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X  function  rbc  -  radiating  bondary  condition 
X  written  by:  Lt.  Hugh  Me  Bride  USNR 
X  date  :  Mar  92 

* 

X  rbc  builds  the  matrix  A(i,k)  which  is  required  by  the  radiation 
%  boundary  condition 

function  v  -  rbc( j , l,pm,dxs,kf ) 

t  variables 

<  j,l  -dimensions  of  the  matrix 
X  pm  -  propagating  modes 
a  =  zeros ( j ) ; 

for  i  =  l: j, 

for  k  =  1:1, 

a(i,k)  =  exp(sqrt(-l) *pm*pi*(i-k) *dxs) ; 
end 


end 


a(:,l)  =  - 5*a  ( :  , 1)  ;  a(:,k)  =  .5*a(:,k); 
v  =  sqrt(kf  ,'2-(pm*pi)  *2)  *a; 
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%  function  rho-  del  funtion  used  for  thr  trapezoidal  rule 
%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

I  generates  a  vector  containing  constants  used  repeatedly  throughout 
%  the  program. 

function  rho  -  rhov(kl,kt,dxs,dts) ; 


rhoi  ~(dts*2)/(kl*2*dxs‘2) ; 
rho2  *=  (dts'2)  /  (kt-'2*dxs~2)  ; 

rho3  =  (dts'2/(4*dxs"2))*((l/kl*2)-(l/kt''2)); 
rho  »=  [rhoi  rho2  rho3],- 
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%  function  stop  -  defines  the  stopping  criteria  for  constructing  the  radiation 
%  boundary  condition  i.e  if  the  result  of  stop  is  0  then  only  the  fundamental  mo 
%  is  allowed  to  propagate,  J  only  the  fundamental  and  the  first  mode  are  allowed 
%  to  prpogate  and  so  on. 

%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Mar  92 

% 

function  v  *  stop(k) 


%  variables 
%  k  - 

%  n  -  nth  propagating  mode 
n  *  0  ;  u  *  -1; 
while,  (k'2  -  n*2*pi~2)  >  0. 
m  *=  m+1;  n  *=  n+1; 

end 
v  =  m; 
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%  function  tcori3  -  treatment  of  corners  l  and  3 
%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

* 

% 

%  tcor!3  applies  the  elastic  wave  equation  at 
%  corners  1  and  3  as  per  fig(  )  since  they  use  the 
%  same  difference  formula. 

%  the  corner  nodes  are  located  (l  first  then  3) 

%  identified  as  p  and  q  (pi  ,ql  in  the  case  of  corner  3) 

%  the  necessary  neighbours  are  picked  off  from  the  u  and  v 
%  matrices  and  weighted  accordingly  and  the  new  updated  values 
%  for  u  and  v  are  computed. 


function  [un,vn]  =  tcorl3 (uc,vc,uo,vo,un, vn, rows, pcoll, pcolr, rho) ; 


%  variables 

%  un  :  updated  values  of  u  vn  :  updated  values  of  v 

%  uc  :  current  values  of  u  vc  :  current  values  of  v 

%  uo  :  old  values  of  u  vo  :  old  values  of  v 

%  rho  :  vector  containing  global  constants 

%  rows  :  used  to  identify  the  elements  of  the  matrix  which 
\  are  zero  for  all  times  they  also  contain  the  row  location  of 
%  the  corner  nodes. 

%  pclor  and  pcoll  carries  out  the  same  function  as  rows  for 
%  the  columns,  and  the  contain  the  column  location  of  the 
%  corner  nodes. 


[il,jl]  *  size(rows); 
ii2,i2j  «  size(pcoir); 

[ i3 , j  3 1  =  size (pcoll) ; 

t  we  pick  off  the  elments  of  rows  and  pcolr  which 
%  identify  the  location  of  corner  1. 

p  *=  rows(il)-l; 
q  =  pcoll (j 3) ; 

*  generate  any  required  local  constants 

cc  «  (2  -  2*rho(l) -2*rho(2) )  ; 

*  tne  updated  values  of  the  corner  nodes  for  u  and  v  are 
%  calculated 


un(p,q)  =  cc*uc(p,q)  -  uo(p,q)... 

+  rho (1) * (uc(p,q-l)  +  uc (p,q+l) ) . . . 

+  rho(2) * (uc (p-1 ,q)  +  uc (p+l , q) ) . . . 

-  2+rho(3)*(vc(p-l,q-l)  +  vc(p+l,q+l)  +  2*vc (p,q) ) . . . 

■'  2*rho  (3)  *  (vc  (p,  q+l)  +  vc(p,q-l)  +vc(p+l,q)  +  vc (p-1, q) ) ; 


vn(p,q)  =  cc*vc(p,q) 

+  rho (2) * (vc (p, q-1)  + 
+  rho(l)*(vc(p-l,q) 


vo(p,q) . . . 
vc(p,q+l) ) . . . 

+  vc (p+l ,q) ) . . . 
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-  2*j.bo(3)  *(uc(p-l,q-l)  +  uctp+ljq-^l)  +  2*uc  (p,q) )  .  . , 

+  2*rho(3)*(uc(p,q+.i)  +  uc(p,q-l)  +  uc(p+l,q)  +  uc(p-l,q)); 


%  the  same  process  is  repeated  for  corner  4  below 


pi  ■  rows(jl)+l; 
ql  -  pcolr(i2) ; 


un(pl,ql)  «cc*uc(pl,ql)  -  uo(pl,ql) _ 

+  rho(l)*(uc(pl,qi-i)  +  uc(pl,ql+l; ) — 

+  rho(2)*(uc(pl-l,ql)  +  uc(pl+l,ql) ) . . . 

-  2*rho(3)*(vc(pl-l,ql-J.)  +  vc(pl+l,ql+l)  +  2*vc(pl,ql) ) . , . 

+  2*rho(3)»(vc(pl,ql+l)  +  vc(pl,ql-l)  +  vc(pi+l,qi)  +  vc(pl-l,ql) ) ; 


vn(pl,ql)  «cc*vc(pl,ql)  -  vo(pl.ql) _ 

+  rho(2) *(vc(plfql-l)  +  vc (pi , ql+1) ) . . . 

+  rho(i)*(vc(pi~i,qi)  +  vc(pl+l,ql) ) . . . 

-  2*rho (3)*{uc(pl-l,ql-l) +uc (pl+1 , ql+1) +2*uc (pl,ql) ) . . . 

+  2*rho(3) * (uc(pl,ql+l)  +  uc(pl,ql-l)  +  uc(pl+i,qi)  +  uc(pl-l,ql) ) ; 
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%  function  tcor24  -  treatment  of  corners  2  and  4 
%  written  by:  Lt.  Hugh  Me  Bride  USNR 

*  date  :  Apr  92 

% 

% 

%  tcor24  applies  the  elastic  wave  equation  at 
%  corners  2  and  4  as  per  fig(  )  since  they  use  the 
%  same  difference  formula. 

%  the  corner  nodes  are  located  (2  first  then  4) 

%  identified  as  p  and  q  (pi  ,ql  in  the  case  of  corner  4) 

%  the  necessary  neighbours  are  picked  off  from  the  u  and  v 
%  matrices  and  weighted  accordingly  and  the  new  updated  values 

*  for  u  and  v  are  computed. 


function  [un,vn]  •  tcor24 (uc,vc,uo,vo, un.vn, rows, pcoll, pcolr, rho) 
*  variables 

%  un  :  updated  values  of  u  vn  :  updated  values  of  v 

%  uc  ;  current  values  of  u  vc  :  current  values  of  v 

%  uo  :  old  values  of  u  vo  :  old  values  of  v 

%  rho  :  vector  containing  global  constants 

%  rows  :  used  to  identify  the  elements  of.  the  matrix  which 
%  are  zero  for  air  times  they  also  contain  the  row  location  cif 
.  the  corner  nodes. 

%  pclor  and  pcoll  carries  out  the  same  function  as  rows  for 
%  the  columns,  and  the  contain  the  column  location  of  the 
%  corner  nodes . 


[il,jl]  *  size ( rows) ; 

[i2,j2]  >-•  •-e(pcolr); 

[  i.3 ,  j 3  2  r pcoll)  ; 

%  we  pick  off  t . ients  of  rows  and  pcolr  which 

%  identify  the  location  of  corner  2. 

p  “  rows ( il) -1; 
q  =  pcolr (i2) ; 

%  generate  any  required  local  constants 
cc  <=  (2  -  2*rho(l)  -2* rho (2) )  ; 


%  the  updated  values  of  the  corner  nodes  for  u  and  v  are 
%  calculated 

un(p,q)  =  cc*uc (p. q)  -  uo(p,q)... 

+  rhO(l)*(uc'p,q-l)  +  uc (p , q+l) ) . . . 

+  rho (2) * (uc (p-1 ,q)  +  uc (p+1 , q> ) . . . 

+  2*rho(3) * (vc(p+l,q-l)  +  vc(p-l,q+l)  +  2 *vc (p , q) ) . . . 

-  2* rho (3) * (vc(p,q+l)  +  vc(p,q-l)  f  vc(p+l,q)  +  vc(p-l,q)); 


vn(p,q)  «=  oc*vc(p,q)  -  vo  ,  q)  .  .  . 

+  rho(2) * (vc (p,q-l)  +  vc (p , q+l ) ) . . . 

+  rho(  1)  *  ( vc  (p-1 ,  q)  +  vc  (pO. ,  q)  )  .  .  . 

+  2*rho(3) *(uc(pvi,q-i)  +  uc(p-l,q+l)  +  2*uc(p,q) ) . . . 

-  2*lho(3) * (uc(p,q+l)  +  uc (p , q-1)  +  uc(p+l,q)  +  uc(p-l,q)); 


%  the  ssine  process  is  repeated  for  corner  4  below 


pi  -  rows(jl)+l; 
ql=  pcoll( j  3 ) ; 


un(pl,ql)  «  cc*uc(pl,ql)  -  uo(pl,ql) — 

+  rho(l) * (uc(pl,ql-l)  +  uc(pl,ql+l) ) . . . 

+  rho(2) *(uc(pl-l,ql)  +  uc(pl+l,ql) ) . . . 

+  2*rhc(3)*(vc(pl+l,ql-l)  +  vc(pl-l,ql+l)  +  2*vc (pl,ql) ) . . . 

-  2*rho(3)*(vc(pl,ql+l)  +  vc (pi , ql-1)  +  vc(pl+l,ql)  +  vc(pl-l,ql) ) ; 


vn(pl,ql)  *  cc*vc(pl,ql)  -  vo(pl,ql)... 

+  rho(2)*(vc(pl,ql-l)  +  vc(pl,ql+l) ) . . . 

+  rho(l)*(vc(pl-l,ql)  +  vc(pl+l,ql) ) — 

+  2*rho(3)*(uc(pl+l,ql-l)  +  uc(pl-l,ql+l)  +  2*uc(pl,ql) ) . . . 

-  2*rho(3) *(uc(pl,ql+l)  +  uc(pl,ql-l)  +  uc(pl+l,ql)  +  uc(pl-l,ql) ) ; 
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%  function  trid  -  tridiagonal 
%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Mar  92 

% 

%  trid  generates  a  tridiagonal  matrix  for  the 
%  radiation  boundary  condition  of  the  fluid  . 

%  the  elements  of  the  sub  and  super  diagonal  are  (dts/kf *dxs) *2 
%  the  main  diagonal  component  is  2  -  4 (dts/kf *dxs) *2 
%  the  (n,2)  and  (l,n-l)  contain  (dts/kf *dxs) *2  to  satisfy  the 
%  periodic  boundary  conditions. 

function  m  «=  trid(dxs,dts,n,kf) 

%  variables 

%  dxs  :  scaled  spacing 

%  dts  i  scaled  time  step 

%  n  :  dimension  of  matrix 

%  kf  :  scaled  constant 


\  an  identity  matrix  for  the  elements  of  the  main  diagonal 
dl  «  eya(n) ; 

t  the  sub  and  super  diagonals 

d2  =  di?:g(cnes(n“  1,  J.) , i )  +  diag(ones(n-l,l)  ,-l) ; 
%  the  required  co-efficients 

lho  *=  dts/  (kf *dxs) ;  rho  «  (rho~2)  ; 

%  generates  the  required  matrix 

d  =2*dl  -4*rho*di  +rho*d2;  t 

d(i,n-l)«=  rho  ;  d(n,2)  =  rho; 

m  ■-  d; 
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%  function  ucnx  -  u  coeff icients/constants  for  the  boundary  facing  the 
%  negative  x  direction. 

%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

%  provides  the  coefficients  developed  by  applying  the  Ilan  and  Lowenthal 
%  technique  to  the  boundary  with  unit  normal (-1,0)  for  the  corresponding  u 
%  and  v  values. 

t 

function  cunx  *  ucnx(kl,kt,dxs,dts) 

%  variables 

%  kl  -  scaled  longitudinal  speed 
%  kt  -  scaled  transverse  speed 
%  dxs  -  scaled  spacing 
%  dts  -  scaled  time  step 

%  the  coefficients  are  calculated  and  stored  in  the  vector  cunx 
%  for  use  in  bndfnx 


cl  rr  2  -  2*(dts*2/dxs*2) *( (l/kl-2) +(l/kt^2) ) ; 
c2  «  2*dts*2/ (dxsA2*kl'2) ; 
c3  =  (dtsA2) / (dxsA2*ktA2)  ; 

c4  «  -  (dts  "'2/  (2*dxsA2 )  )*(  { l/kl*2) - (  l/ktA2)  )  ; 
c5  «  (dtsA2/ (2*dxsA2) ) * ( (l/klA2)-(3/ktA2) ) ; 
cunx  =  [cl  c2  C3  c4  c5 ] ; 
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%  function  ucny  -  u  coefficients/constants  for  the  boundary  facing  the 
%  negative  y  direction. 

*  written  by:  Lt.  Hugh  Me  Bride  USNE 
%  date  :  Apr  9? 

i  provides  the  coefficients  developed  by  applying  the  llan  and  Lowenthal 
?■  technique  to  the  boundary  with  unit  normal  (0, -1)  for  the  corresponding  u 
x  and  v  values, 
t 


function  cuny  *  ucny(kl,kt,dxs,dts) 


*  vai iables 

%  kl  -  scaled  longitudinal  speed 
%  kt  ~  scaled  transverse  speed 
%  dxs  -  scaled  spacing 
V  dts  -  scaled  time  step 

%  the  coefficients  are  calculated  and  stored  in  the  vector  cuny 
%  for  use  in  bndfny 

cl  *  2  -  2*(dts*2/dxs*2)*((l/kl''2)+(l/kt'2))  ; 

C2  *  2*(dts'2/ (dxs'2*kt'2) )  ; 

C3  *  (dts'2)  /  (dxs'2*kl'2)  ; 

c4  «  - (dtr.'2/  (2*dxsA2) )  *(  (i/kl'2)-(l/kt'2) )  ; 
c5  «  (dts' 2/  (2*dxs“2)  )*((3/kt*'2)-(l/klA2))  ; 
cuny  *  [cl  c2  c3  c4  cS] ; 
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%  function  ucpy  ■  u  coef f icients/constants  for  the  boundary  facing  the 
%  positive  x  6.  ."sot  ion. 

%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

%  provides  the  coefficients  developed  by  applying  the  Ilan  and  Lowenthal 

%  technique  to  the  boundary  with  unit  normal (1,0)  for  the  corresponding  u 

%  and  v  values. 

% 

function  cupx  -  ucpx(kl,kt,dxs,dts) 

%  variables 

%  kl  -  sc.altd  longitudinal  epeed 
%  kt  -  scaled  transverse  speed 
%  dxs  -  scaled  spacing 

%  dts  -  scaled  tine  step 

%  the  coefficients  are  calculated  and  stored  in  the  vector  cupx 
%  for  use  in  bndfpx 


Cl  =  2  -  2*  (dts "2 /dxs '2)  *  (  ( i/kl"2) +(l/kt"2) ) ; 
c2  ■=  2*(dts"2/(dxsA2*kl"2)  )  ; 
c3  ■=  (dts"2) / (dxs"2*kt"2) ; 
c4  =  (dts" 2/  (2*dxs"2)  )*(  (i/kl"2) -(l/kt*2) )  ; 
c5  =  - (dts "2/ (2*dxs"2) )*((l/kl"2)-(3/kt"2)) ; 
cupx  =  [cl  c2  c3  c4  c5] ; 
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%  function  ucny  -  u  coefficients/constants  for  the  boundary  facing  the 
l  positive  y  direction. 

%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

%  provides  the  coefficients  developed  by  applying  the  Ilan  and  Lowenthal 

%  technique  to  the  boundary  with  unit  normal (0,1)  for  the  corresponding  u 

%  and  v  values. 

% 

function  cupy  -  uepy (kl,kt,dxs,dts) 

%  variables 

%  kl  -  scaled  longitudinal  speed 
%  kt  -  scaled  transverse  speed 
%  dxs  -  scaled  spacing 

%  dts  -  scaled  time  step 

%  the  coefficients  are  calculated  and  stored  in  the  vector  cupy 
%  for  use  in  bndfpy 


cl  =  2  -  2*(dts"2/dxs-'2)*(  (l/kl'2) +(l/kt'2)  )  ; 
c2  »  2*(dts"2/ (dxs'2*kt'2) ) ; 
c3  =  (dts *2) / (dxsA2*kl*2) ; 

04  =  (dts~2/(2*dxs~2))*((l/k)'2)-(l/kt-2)) ; 
cj  »  (dts‘2/(2*dxs*2))*((l/kl'‘2)-(3/kf'2))  ; 
cupy  ■*  (cl  C2  c3  C4  c5]  ; 
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i  function  vcnx  -  v  coef'icients/coristants  for  the  boundary  facing  the 
%  negative  x  direction. 

%  written  by:  Lt.  Hugh  fc.c  Bride  USNR 
\  date  :  Apr  92 

%  provides  the  coefficients.  developed  hy  applying  the  Ilan  and  Loventhal 
%  technique  to  the  boundary  with  unit  normal (-1,0)  for  the  corresponding  u 
%  and  v  values. 

t 

function  cvnx  -  vcnx(kl,kt,dxs,dts) 


%  variables 

%  kl  -  scaled  longitudinal  speed 
%  kt  -  scaled  transverse  speed 
%  dxs  -  scaled  spacing 
%  dts  -  scaled  time  step 

%  the  coefficients  are  calculated  and  stored  in  the  vector  cvnx 
4  for  use  in  bndfnx 


dl  =  2  -  2*(dts-'2/dxs"2)  *  ( (l/kl*2) +(l/kt*2)  )  ; 
d2  -  2*(dtsA2/(dxs  '2*ktA2))  ; 
d3  =  (dts'2) / (dXSA2*klA2) ; 

d4  =  -<dtsA2/(2*dxsA2) )*( (l/klA2)-(l/ktA2) ) ; 
d5  «  (dtsA2/(2*dxs^2))*((3/kt*2)-(l/kl-2)); 
cvnx  “  [dl  d2  d3  d4  d5] ; 
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%  function  vcny  -  v  coefficients/constants  for  the  boundary  facing  the 
%  negative  y  direction. 

%  written  by:  Lt.  Hugh  He  Bride  USNR 
%  date  :  Apr  92 

%  provides  the  coefficients  developed  by  applying  the  Ilan  and  Lowenthal 
%  technique  to  the  boundary  with  unit  normal (o, -1)  for  the  corresponding  u 
%  and  v  values. 

* 


function  cvny  ■  vcny (kl,kt,dxs,dts)  * 

\  variables 

%  kl  -  scaled  longitudinal  speed 
%  kt  -  scaled  transverse  speed 
%  dxs  -  scaled  spacing 
$  dts  -  scaled  time  step 

%  the  coefficients  are  calculated  and  stored  in  the  vector  cvny 
%  for  use  in  bndfny 


dl  =  2  -  2*  (dts *2 /dxsA2)  *((l/klA2)  +  (l/kt'>2))  ; 

d2  =  2*(dts',2/(dxs''2*ki'2))  J 

d3  =  (dts *2 )  /  (dxs',2*kt*2)  ; 

d4  =  -(dts*2/ (2*dxs-2) )*( ( l/kl*2) - ( l/kt'2 ) ) ; 
dS  *  (dts“2/(2*dxs*2) )*( (l/kl'2)-(3/kt*2) ) ; 
cvny  =  [dl  d2  d3  dl  d5J; 


A 
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%  function  vcpx  -  v  coefficients/constants  for  the  boundary  facing  the 
%  positive  x  direction. 

%  written  by:  Lt.  Hugh  Me  Bride  USNR 
%  date  :  Apr  92 

%  provides  the  coefficients  developed  by  applying  the  Ilan  and  Lowenthal 
%  technique  to  the  boundary  with  unit  normal(i,0)  for  the  corresponding  u 
%  and  v  values. 

% 

function  cvpx  -  vcpx(kl,kt,dxs,dtp) 

%  variables 

%  kl  -  scaled  longitudinal  speed 
%  kt  -  scaled  transverse  speed 
%  dxs  -  scaled  spacing 
%  dts  -  scaled  time  step 

%  the  coefficients  are  calculated  and  stored  in  the  vector  cvpx 
%  for  use  in  bndfpx 


dl  -  2  -2* (dts*2/dxs'2) *( (l/klA2)+(l/kt*2) )  ; 
d2  ■*  2 *  ( dts “2/  { dxs ' 2 *kt ~ 2 ) ) ; 
d3  *=  (dts'2)  /  (dxs  '2*kl*2)  ; 

d4  >=  (dts*2/(2*dxs*2))*((l/kl*2)-(l/kt*2)); 
d5  =  (dts'2/(2*dxs''2))*((i/kl,'2)-(3/kr.-2)); 
cvpx  ■>  [dl  d2  d3  d4  d5]  ; 


%  function  vcny  -  v  coefficients/constants  for  the  boundary  facing  the 
1  negative  y  direction. 

X  written  by:  Lt.  Hugh  He  Bride  USNR 
%  data  :  Apr  92 

*  provides  the  coefficients  developed  by  applying  the  Ilan  and  Lowenthal 
X  technique  to  the  boundary  with  unit  norma  1(0,1)  for  the  corresponding  u 
X  and  v  values. 

X 


function  cvpy  =  vcpy(k.l,kt,dxs,dts)  t 

%  variables 

*  kl  •-  scaled  longitudinal  speed 
X  kt  -  scaled  transverse  speed 

%  dxs  -  scaled  spacing 
X  dts  -  scaled  time  step 

*  the  coefficients  are  calculated  and  stored  in  the  vector  cvpy 
%  for  use  in  bndfpy 


dl  *  2  -  2* (dts"2/dxs"2 ) * ( (i/kl*2) + (l/kt'2) ) ; 
d2  =  2* (dts"2/ (dxs"2*kl"2) )  ; 
d3  »  (dts"2) / (dxs*2*kt*2)  ; 
d4  =  (dts"  2  /  (2*dxs"2))*((l/kl"2)-(l/kt"2))  ; 
d5  «  -(dts"2/(2*dxs"2) ) *  ( (l/kl"2)-(3/kt"2) )  ; 
cvpy  =  [dl  d2  d3  d4  d5]; 


■6 


A 


126 


REFERENCES 


Fuyuki,  M.,  and  Matsumoto,  Y.,  "Finite  Difference  Analysis  of  Rayleigh 
Wave  Scattering  at  a  Trench/'  Bulletin  of  the  Seismological  Society  of 
America,  v.  70,  December  1980. 

Ilan  A.,  and  Lowenthal,  D.,  "instability  of  Finite  Difference  Schemes  due 
to  Boundary  Conditions  in  Elastic  Media,"  Geophysical  Prospecting,  v.  24, 
1976. 

Interview  between  Prof.  C.  L.  Scandrett,  Department  of  Mathematics, 
Naval  Postgraduate  School,  and  the  author,  30  May,  1992. 

Kinsler,  L.  E.,  Frey,  A.  R.,  Coppens,  A.  B.,  Sanders,  J.  V.,  Fundamentals  of 
Acoustics,  Third  Edition,  John  Wiley  &  Sons,  1982. 

Lalanne  M.,  Berther  P.,  and  Der  Hagopian,  J.,  Mechanical  Vibration  for 
Engineers,  First  Edition, John  Wiley  &  Sons,  1982. 

Scandrett,  C.  L.,  and  Kreigsmann,  G.  A.,  "Decoupling  Approximation 
Applied  to  an  Infinite  Array  of  Fluid  Loaded  Baffled  Membranes,"  to  be 
submitted  for  publication. 


127 


INITIAL  DISTRIBUTION  LIST 


* 


1 .  Library  (Code  52) . 2 

Naval  Postgraduate  School 

Monterey,  CA  93943-5000 

2.  Defense  Technical  Information  Center . 2 

Cameron  Station 

Alexandria,  VA  22304-6145 

3.  Prof.  Clyde  Scandrett . 1 

Code  MA/Sd 

Department  of  Mathematics 
Naval  Postgraduate  School, 

Monterey,  CA  93943-5000 

4.  Prof.  Van  Emden  Henson . 1 

Code  MA/Hv 

Department  of  Mathematics 
Naval  Postgraduate  School, 

Monterey,  CA  93943-5000 

5.  Prof.  Young  W.  Kwon . 1 

Code  ME/Kw 

Department  of  Mechanical  Engineering 
Naval  Postgraduate  School, 

Monterey,  CA  93943-5000 


6.  LT  Hugh  McBride . 2 

201  Glenwood  Circle,  Apt  26D 

Monterey,  CA  93940 

7.  Chairman . 1 

Department  of  Mathematics 

Naval  Postgraduate  School 
Monterey,  CA  93943 


4 


128 


